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 9490 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2018-04-23T10:44:07+02:00 (6 years ago)
Author:
gm
Message:

#2075 - dev_merge_2017: scale-aware setting of lateral viscous and diffusive coefficient

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90

    r9125 r9490  
    5353         WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics' 
    5454         WRITE(numout,*) '~~~~~~~~~~~~' 
    55          WRITE(numout,*) 'Namelist nam_diatmb : set tmb outputs ' 
    56          WRITE(numout,*) 'Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb 
     55         WRITE(numout,*) '   Namelist nam_diatmb : set tmb outputs ' 
     56         WRITE(numout,*) '      Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb 
    5757      ENDIF 
    5858      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r9436 r9490  
    148148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
    149149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
    150  
    151150 
    152151   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r9405 r9490  
    316316      IF(lwp) THEN                  ! control print 
    317317         WRITE(numout,*) '   Namelist : namrun   ---   run parameters' 
    318          WRITE(numout,*) '      job number                      nn_no           = ', nn_no 
     318         WRITE(numout,*) '      Assimilation cycle              nn_no           = ', nn_no 
    319319         WRITE(numout,*) '      experiment name for output      cn_exp          = ', TRIM( cn_exp           ) 
    320320         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in    = ', TRIM( cn_ocerst_in     ) 
     
    351351      ENDIF 
    352352 
    353       no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon) 
    354       cexper = cn_exp 
     353      cexper = cn_exp         ! conversion DOCTOR names into model names (this should disappear soon) 
    355354      nrstdt = nn_rstctl 
    356355      nit000 = nn_it000 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r9168 r9490  
    99   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module  
    1010   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread 
    11    !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys 
     11   !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + remove CPP keys 
    1212   !!---------------------------------------------------------------------- 
    1313 
     
    2929   PUBLIC   dta_tsd        ! called by istate.F90 and tradmp.90 
    3030 
    31    LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag 
    32    LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag 
     31   !                                  !!* namtsd  namelist : Temperature & Salinity Data * 
     32   LOGICAL , PUBLIC ::   ln_tsd_init   !: T & S data flag 
     33   LOGICAL , PUBLIC ::   ln_tsd_dmp    !: internal damping toward input data flag 
    3334 
    3435   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     
    5859      TYPE(FLD_N)                   ::   sn_tem, sn_sal 
    5960      !! 
    60       NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal 
     61      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal 
    6162      !!---------------------------------------------------------------------- 
    6263      ! 
     
    7273      IF(lwm) WRITE ( numond, namtsd ) 
    7374 
    74       IF( PRESENT( ld_tradmp ) )   ln_tsd_tradmp = .TRUE.     ! forces the initialization when tradmp is used 
     75      IF( PRESENT( ld_tradmp ) )   ln_tsd_dmp = .TRUE.     ! forces the initialization when tradmp is used 
    7576       
    7677      IF(lwp) THEN                  ! control print 
     
    7980         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    8081         WRITE(numout,*) '   Namelist namtsd' 
    81          WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init 
    82          WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp 
     82         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init = ', ln_tsd_init 
     83         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp  = ', ln_tsd_dmp 
    8384         WRITE(numout,*) 
    84          IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_tradmp ) THEN 
     85         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN 
    8586            WRITE(numout,*) 
    86             WRITE(numout,*) '   T & S data not used' 
     87            WRITE(numout,*) '   ===>>   T & S data not used' 
    8788         ENDIF 
    8889      ENDIF 
     
    9596      ! 
    9697      !                             ! allocate the arrays (if necessary) 
    97       IF(  ln_tsd_init .OR. ln_tsd_tradmp ) THEN 
     98      IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 
    9899         ! 
    99100         ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     
    129130      !!              - 'key_orca_lev10' interpolates on 10 times more levels 
    130131      !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh 
    131       !!              - ln_tsd_tradmp=F: deallocates the T-S data structure 
     132      !!              - ln_tsd_dmp=F: deallocates the T-S data structure 
    132133      !!                as T-S data are no are used 
    133134      !! 
     
    149150      ! 
    150151      !                                   !==   ORCA_R2 configuration and T & S damping   ==!  
    151       IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_tradmp ) THEN    ! some hand made alterations 
     152      IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_dmp ) THEN    ! some hand made alterations 
    152153         ! 
    153154         ij0 = 101   ;   ij1 = 109                       ! Reduced T & S in the Alboran Sea 
     
    238239      ENDIF 
    239240      ! 
    240       IF( .NOT.ln_tsd_tradmp ) THEN                   !==   deallocate T & S structure   ==!  
     241      IF( .NOT.ln_tsd_dmp ) THEN                   !==   deallocate T & S structure   ==!  
    241242         !                                              (data used only for initialisation) 
    242243         IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run' 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r9190 r9490  
    9898      !!---------------------------------------------------------------------- 
    9999      ! 
     100      IF(lwp) THEN 
     101         WRITE(numout,*) 
     102         WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme' 
     103         WRITE(numout,*) '~~~~~~~~~~~~' 
     104      ENDIF 
     105      ! 
    100106      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 
    101107      READ  ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) 
     
    107113 
    108114      IF(lwp) THEN                    ! Namelist print 
    109          WRITE(numout,*) 
    110          WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme' 
    111          WRITE(numout,*) '~~~~~~~~~~~~' 
    112115         WRITE(numout,*) '   Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 
    113116         WRITE(numout,*) '      linear dynamics : no momentum advection          ln_dynadv_NONE = ', ln_dynadv_NONE 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r9168 r9490  
    5353   PUBLIC   dyn_hpg_init   ! routine called by opa module 
    5454 
    55    !                                 !!* Namelist namdyn_hpg : hydrostatic pressure gradient 
    56    LOGICAL , PUBLIC ::   ln_hpg_zco   !: z-coordinate - full steps 
    57    LOGICAL , PUBLIC ::   ln_hpg_zps   !: z-coordinate - partial steps (interpolation) 
    58    LOGICAL , PUBLIC ::   ln_hpg_sco   !: s-coordinate (standard jacobian formulation) 
    59    LOGICAL , PUBLIC ::   ln_hpg_djc   !: s-coordinate (Density Jacobian with Cubic polynomial) 
    60    LOGICAL , PUBLIC ::   ln_hpg_prj   !: s-coordinate (Pressure Jacobian scheme) 
    61    LOGICAL , PUBLIC ::   ln_hpg_isf   !: s-coordinate similar to sco modify for isf 
    62  
    63    INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     55   !                                !!* Namelist namdyn_hpg : hydrostatic pressure gradient 
     56   LOGICAL, PUBLIC ::   ln_hpg_zco   !: z-coordinate - full steps 
     57   LOGICAL, PUBLIC ::   ln_hpg_zps   !: z-coordinate - partial steps (interpolation) 
     58   LOGICAL, PUBLIC ::   ln_hpg_sco   !: s-coordinate (standard jacobian formulation) 
     59   LOGICAL, PUBLIC ::   ln_hpg_djc   !: s-coordinate (Density Jacobian with Cubic polynomial) 
     60   LOGICAL, PUBLIC ::   ln_hpg_prj   !: s-coordinate (Pressure Jacobian scheme) 
     61   LOGICAL, PUBLIC ::   ln_hpg_isf   !: s-coordinate similar to sco modify for isf 
     62 
     63   !                                !! Flag to control the type of hydrostatic pressure gradient 
     64   INTEGER, PARAMETER ::   np_ERROR  =-10   ! error in specification of lateral diffusion 
     65   INTEGER, PARAMETER ::   np_zco    =  0   ! z-coordinate - full steps 
     66   INTEGER, PARAMETER ::   np_zps    =  1   ! z-coordinate - partial steps (interpolation) 
     67   INTEGER, PARAMETER ::   np_sco    =  2   ! s-coordinate (standard jacobian formulation) 
     68   INTEGER, PARAMETER ::   np_djc    =  3   ! s-coordinate (Density Jacobian with Cubic polynomial) 
     69   INTEGER, PARAMETER ::   np_prj    =  4   ! s-coordinate (Pressure Jacobian scheme) 
     70   INTEGER, PARAMETER ::   np_isf    =  5   ! s-coordinate similar to sco modify for isf 
     71   ! 
     72   INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
    6473 
    6574   !! * Substitutions 
     
    95104      ! 
    96105      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    97       CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    98       CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
    99       CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
    100       CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    101       CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    102       CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
     106      CASE ( np_zco )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
     107      CASE ( np_zps )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
     108      CASE ( np_sco )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
     109      CASE ( np_djc )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
     110      CASE ( np_prj )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
     111      CASE ( np_isf )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
    103112      END SELECT 
    104113      ! 
     
    179188      ENDIF 
    180189      ! 
    181       !                               ! Set nhpg from ln_hpg_... flags 
    182       IF( ln_hpg_zco )   nhpg = 0 
    183       IF( ln_hpg_zps )   nhpg = 1 
    184       IF( ln_hpg_sco )   nhpg = 2 
    185       IF( ln_hpg_djc )   nhpg = 3 
    186       IF( ln_hpg_prj )   nhpg = 4 
    187       IF( ln_hpg_isf )   nhpg = 5 
    188       ! 
    189       !                               ! Consistency check 
     190      !                               ! Set nhpg from ln_hpg_... flags & consistency check 
     191      nhpg   = np_ERROR 
    190192      ioptio = 0 
    191       IF( ln_hpg_zco )   ioptio = ioptio + 1 
    192       IF( ln_hpg_zps )   ioptio = ioptio + 1 
    193       IF( ln_hpg_sco )   ioptio = ioptio + 1 
    194       IF( ln_hpg_djc )   ioptio = ioptio + 1 
    195       IF( ln_hpg_prj )   ioptio = ioptio + 1 
    196       IF( ln_hpg_isf )   ioptio = ioptio + 1 
     193      IF( ln_hpg_zco ) THEN   ;   nhpg = np_zco   ;   ioptio = ioptio +1   ;   ENDIF 
     194      IF( ln_hpg_zps ) THEN   ;   nhpg = np_zps   ;   ioptio = ioptio +1   ;   ENDIF 
     195      IF( ln_hpg_sco ) THEN   ;   nhpg = np_sco   ;   ioptio = ioptio +1   ;   ENDIF 
     196      IF( ln_hpg_djc ) THEN   ;   nhpg = np_djc   ;   ioptio = ioptio +1   ;   ENDIF 
     197      IF( ln_hpg_prj ) THEN   ;   nhpg = np_prj   ;   ioptio = ioptio +1   ;   ENDIF 
     198      IF( ln_hpg_isf ) THEN   ;   nhpg = np_isf   ;   ioptio = ioptio +1   ;   ENDIF 
     199      ! 
    197200      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    198201      !  
     202      IF(lwp) THEN 
     203         WRITE(numout,*) 
     204         SELECT CASE( nhpg ) 
     205         CASE( np_zco )   ;   WRITE(numout,*) '   ==>>>   z-coord. - full steps ' 
     206         CASE( np_zps )   ;   WRITE(numout,*) '   ==>>>   z-coord. - partial steps (interpolation)' 
     207         CASE( np_sco )   ;   WRITE(numout,*) '   ==>>>   s-coord. (standard jacobian formulation)' 
     208         CASE( np_djc )   ;   WRITE(numout,*) '   ==>>>   s-coord. (Density Jacobian: Cubic polynomial)' 
     209         CASE( np_prj )   ;   WRITE(numout,*) '   ==>>>   s-coord. (Pressure Jacobian: Cubic polynomial)' 
     210         CASE( np_isf )   ;   WRITE(numout,*) '   ==>>>   s-coord. (standard jacobian formulation) for isf' 
     211         END SELECT 
     212         WRITE(numout,*) 
     213      ENDIF 
    199214      !                           
    200215      IF ( .NOT. ln_isfcav ) THEN     !--- no ice shelf load 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r9190 r9490  
    1717   USE phycst         ! physical constants 
    1818   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    19    USE ldfslp         ! lateral diffusion: slopes of mixing orientation 
    2019   USE dynldf_lap_blp ! lateral mixing   (dyn_ldf_lap & dyn_ldf_blp routines) 
    2120   USE dynldf_iso     ! lateral mixing                 (dyn_ldf_iso routine ) 
     
    3433   PUBLIC   dyn_ldf       ! called by step module  
    3534   PUBLIC   dyn_ldf_init  ! called by opa  module  
    36  
    37    !                      ! Parameter to control the type of lateral viscous operator 
    38    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   !: error in setting the operator 
    39    INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   !: without operator (i.e. no lateral viscous trend) 
    40    !                          !!      laplacian     !    bilaplacian    ! 
    41    INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator 
    42    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator 
    43  
    44    INTEGER, PUBLIC ::   nldf   !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
    4535 
    4636   !! * Substitutions 
     
    7262      ENDIF 
    7363 
    74       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
     64      SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend 
    7565      ! 
    76       CASE ( np_lap   )    ;   CALL dyn_ldf_lap  ( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
    77       CASE ( np_lap_i )    ;   CALL dyn_ldf_iso  ( kt )                         ! rotated      laplacian 
    78       CASE ( np_blp   )    ;   CALL dyn_ldf_blp  ( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
     66      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
     67      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian 
     68      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
    7969      ! 
    8070      END SELECT 
     
    10191      !! ** Purpose :   initializations of the horizontal ocean dynamics physics 
    10292      !!---------------------------------------------------------------------- 
    103       INTEGER ::   ioptio, ierr   ! temporary integers  
    104       !!---------------------------------------------------------------------- 
    10593      ! 
    106       !                                !==  Namelist nam_dynldf  ==!   already read in ldfdyn module 
    107       ! 
    108       IF(lwp) THEN                     !== Namelist print  ==! 
     94      IF(lwp) THEN                     !==  Namelist print  ==! 
    10995         WRITE(numout,*) 
    11096         WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics' 
    11197         WRITE(numout,*) '~~~~~~~~~~~~' 
    112          WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 
    113          WRITE(numout,*) '      Type of operator' 
    114          WRITE(numout,*) '         no explicit diffusion       ln_dynldf_NONE = ', ln_dynldf_NONE 
    115          WRITE(numout,*) '         laplacian operator          ln_dynldf_lap  = ', ln_dynldf_lap 
    116          WRITE(numout,*) '         bilaplacian operator        ln_dynldf_blp  = ', ln_dynldf_blp 
    117          WRITE(numout,*) '      Direction of action' 
    118          WRITE(numout,*) '         iso-level                   ln_dynldf_lev  = ', ln_dynldf_lev 
    119          WRITE(numout,*) '         horizontal (geopotential)   ln_dynldf_hor  = ', ln_dynldf_hor 
    120          WRITE(numout,*) '         iso-neutral                 ln_dynldf_iso  = ', ln_dynldf_iso 
    121       ENDIF 
    122       !                                !==  use of lateral operator or not  ==! 
    123       nldf = np_ERROR 
    124       ioptio = 0 
    125       IF( ln_dynldf_NONE ) THEN   ;   nldf = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF 
    126       IF( ln_dynldf_lap  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF 
    127       IF( ln_dynldf_blp  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF 
    128       IF( ioptio /= 1    )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
    129       ! 
    130       IF(.NOT.ln_dynldf_NONE ) THEN    !==  direction ==>> type of operator  ==! 
    131          ioptio = 0 
    132          IF( ln_dynldf_lev )   ioptio = ioptio + 1 
    133          IF( ln_dynldf_hor )   ioptio = ioptio + 1 
    134          IF( ln_dynldf_iso )   ioptio = ioptio + 1 
    135          IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 
     98         WRITE(numout,*) '   Namelist namdyn_ldf: already read in ldfdyn module' 
     99         WRITE(numout,*) '      see ldf_dyn_init report for lateral mixing parameters' 
     100         WRITE(numout,*) 
    136101         ! 
    137          !                             ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 
    138          ierr = 0 
    139          IF( ln_dynldf_lap ) THEN         ! laplacian operator 
    140             IF( ln_zco ) THEN                ! z-coordinate 
    141                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
    142                IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
    143                IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
    144             ENDIF 
    145             IF( ln_zps ) THEN                ! z-coordinate with partial step 
    146                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level              (no rotation) 
    147                IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level              (no rotation) 
    148                IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
    149             ENDIF 
    150             IF( ln_sco ) THEN                ! s-coordinate 
    151                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
    152                IF ( ln_dynldf_hor )   nldf = np_lap_i   ! horizontal             (   rotation) 
    153                IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
    154             ENDIF 
    155          ENDIF 
    156          ! 
    157          IF( ln_dynldf_blp ) THEN         ! bilaplacian operator 
    158             IF( ln_zco ) THEN                ! z-coordinate 
    159                IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
    160                IF( ln_dynldf_hor )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
    161                IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
    162             ENDIF 
    163             IF( ln_zps ) THEN                ! z-coordinate with partial step 
    164                IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation) 
    165                IF( ln_dynldf_hor )   nldf = np_blp     ! iso-level              (no rotation) 
    166                IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
    167             ENDIF 
    168             IF( ln_sco ) THEN                ! s-coordinate 
    169                IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation) 
    170                IF( ln_dynldf_hor )   ierr = 2          ! horizontal             (   rotation) 
    171                IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
    172             ENDIF 
    173          ENDIF 
    174          ! 
    175          IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 
    176          ! 
    177          IF( nldf == np_lap_i )   l_ldfslp = .TRUE.      ! rotation require the computation of the slopes 
    178          ! 
    179       ENDIF 
    180  
    181       IF(lwp) THEN 
    182          WRITE(numout,*) 
    183          IF( nldf == np_no_ldf )   WRITE(numout,*) '   ==>>>   NO lateral viscosity' 
    184          IF( nldf == np_lap    )   WRITE(numout,*) '   ==>>>   iso-level laplacian operator' 
    185          IF( nldf == np_lap_i  )   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background' 
    186          IF( nldf == np_blp    )   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator' 
     102         SELECT CASE( nldf_dyn )             ! print the choice of operator 
     103         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity' 
     104         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator' 
     105         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background' 
     106         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator' 
     107         END SELECT 
    187108      ENDIF 
    188109      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r9124 r9490  
    108108      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars 
    109109      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      - 
    110       REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4            !   -      - 
     110      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      - 
    111111      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace 
    112112      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      - 
     
    139139         ! 
    140140       ENDIF 
    141  
     141          
     142      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max 
     143       
    142144      !                                                ! =============== 
    143145      DO jk = 1, jpkm1                                 ! Horizontal slab 
     
    174176                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
    175177 
    176                   zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     178                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    177179    
    178180                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    & 
     
    189191                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp ) 
    190192 
    191                   zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     193                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    192194 
    193195                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
     
    206208                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
    207209 
    208                zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     210               zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    209211 
    210212               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
     
    227229                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
    228230 
    229                zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     231               zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    230232 
    231233               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    & 
     
    244246                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
    245247 
    246                   zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     248                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    247249 
    248250                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    & 
     
    259261                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    260262 
    261                   zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     263                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    262264 
    263265                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     
    335337         DO jk = 2, jpkm1 
    336338            DO ji = 2, jpim1 
    337                zcof0 = 0.5_wp * rn_aht_0 * umask(ji,jj,jk) 
     339               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk) 
    338340               ! 
    339341               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
     
    353355                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  ) 
    354356               ! vertical mixing coefficient (akzu) 
    355                ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
    356                akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 
     357               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0 
     358               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0 
    357359            END DO 
    358360         END DO 
     
    361363         DO jk = 2, jpkm1 
    362364            DO ji = 2, jpim1 
    363                zcof0 = 0.5_wp * rn_aht_0 * vmask(ji,jj,jk) 
     365               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk) 
    364366               ! 
    365367               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
     
    379381                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  ) 
    380382               ! vertical mixing coefficient (akzv) 
    381                ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
    382                akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 
     383               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0 
     384               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0 
    383385            END DO 
    384386         END DO 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r9250 r9490  
    1919   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
    2020   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form 
    21    USE dynldf    ,ONLY: nldf, np_lap_i   ! dynamics: type of lateral mixing  
    2221   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing  
    23    USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
     22   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator 
    2423   USE trd_oce        ! trends: ocean variables 
    2524   USE trddyn         ! trend manager: dynamics 
     
    156155      !                    !* Matrix construction 
    157156      zdt = r2dt * 0.5 
    158       IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu) 
     157      SELECT CASE( nldf_dyn ) 
     158      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    159159         DO jk = 1, jpkm1 
    160160            DO jj = 2, jpjm1  
     
    171171            END DO 
    172172         END DO 
    173       ELSE                          ! standard case 
     173      CASE DEFAULT               ! iso-level lateral mixing 
    174174         DO jk = 1, jpkm1 
    175175            DO jj = 2, jpjm1  
     
    184184            END DO 
    185185         END DO 
    186       ENDIF 
     186      END SELECT 
    187187      ! 
    188188      DO jj = 2, jpjm1     !* Surface boundary conditions 
     
    274274      !                       !* Matrix construction 
    275275      zdt = r2dt * 0.5 
    276       IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu) 
     276      SELECT CASE( nldf_dyn ) 
     277      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    277278         DO jk = 1, jpkm1 
    278279            DO jj = 2, jpjm1    
    279280               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280281                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    281                   zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     282                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    282283                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    283                   zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     284                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    284285                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    285286                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     
    289290            END DO 
    290291         END DO 
    291       ELSE                          ! standard case 
     292      CASE DEFAULT               ! iso-level lateral mixing 
    292293         DO jk = 1, jpkm1 
    293294            DO jj = 2, jpjm1    
    294295               DO ji = fs_2, fs_jpim1   ! vector opt. 
    295296                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    296                   zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    297                   zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     297                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     298                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    298299                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    299300                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     
    302303            END DO 
    303304         END DO 
    304       ENDIF 
     305      END SELECT 
    305306      ! 
    306307      DO jj = 2, jpjm1        !* Surface boundary conditions 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r9367 r9490  
    2828   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
    2929   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    30    INTEGER       ::   nn_no            !: job number 
    3130   INTEGER       ::   nn_rstctl        !: control of the time step (0, 1 or 2) 
    3231   INTEGER       ::   nn_rstssh   = 0  !: hand made initilization of ssh or not (1/0) 
     
    4645   LOGICAL       ::   ln_xios_read     !: use xios to read single file restart 
    4746   INTEGER       ::   nn_wxios         !: write resart using xios 0 - no, 1 - single, 2 - multiple file output 
     47   INTEGER       ::   nn_no            !: Assimilation cycle 
    4848 
    4949#if defined key_netcdf4 
     
    7474 
    7575   CHARACTER(lc) ::   cexper                      !: experiment name used for output filename 
    76    INTEGER       ::   no                          !: job number 
    7776   INTEGER       ::   nrstdt                      !: control of the time step (0, 1 or 2) 
    7877   INTEGER       ::   nit000                      !: index of the first time step 
     
    122121   INTEGER ::   numstp          =   -1      !: logical unit for time step 
    123122   INTEGER ::   numtime         =   -1      !: logical unit for timing 
    124    INTEGER ::   numout          =    6      !: logical unit for output print; Set to stdout to ensure any early 
    125                                             ! output can be collected; do not change 
     123   INTEGER ::   numout          =    6      !: logical unit for output print; Set to stdout to ensure any 
     124      !                                     !  early output can be collected; do not change 
    126125   INTEGER ::   numnam_ref      =   -1      !: logical unit for reference namelist 
    127126   INTEGER ::   numnam_cfg      =   -1      !: logical unit for configuration specific namelist 
     
    159158 
    160159   !!---------------------------------------------------------------------- 
    161    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     160   !! NEMO/OPA 4.0 , NEMO Consortium (2018) 
    162161   !! $Id$ 
    163162   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90

    r9094 r9490  
    2626   PUBLIC   ldf_c2d   ! called by ldftra and ldfdyn modules 
    2727 
     28   REAL(wp) ::   r1_2  = 0.5_wp           ! =1/2 
     29   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4 
     30   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12 
    2831  
    29    !! * Substitutions 
    30 #  include "vectopt_loop_substitute.h90" 
    3132   !!---------------------------------------------------------------------- 
    3233   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     
    3637CONTAINS 
    3738 
    38    SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 ) 
     39   SUBROUTINE ldf_c1d( cd_type, pahs1, pahs2, pah1, pah2 ) 
    3940      !!---------------------------------------------------------------------- 
    4041      !!                  ***  ROUTINE ldf_c1d  *** 
     
    4344      !! 
    4445      !! ** Method  :   1D eddy diffusivity coefficients F( depth ) 
    45       !!                Reduction by prat from surface to bottom  
     46      !!                Reduction by zratio from surface to bottom  
    4647      !!                hyperbolic tangent profile with inflection point  
    4748      !!                at zh=500m and a width of zw=200m 
     
    5051      !!             DYN      pah1, pah2 defined at T- and F-points 
    5152      !!---------------------------------------------------------------------- 
    52       CHARACTER(len=2)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers 
    53       REAL(wp)                        , INTENT(in   ) ::   prat           ! ratio surface/deep values           [-] 
     53      CHARACTER(len=3)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers 
    5454      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pahs1, pahs2   ! surface value of eddy coefficient   [m2/s] 
    5555      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pah1 , pah2    ! eddy coefficient                    [m2/s] 
     
    5858      REAL(wp) ::   zh, zc, zdep1   ! local scalars 
    5959      REAL(wp) ::   zw    , zdep2   !   -      - 
     60      REAL(wp) ::   zratio          !   -      - 
    6061      !!---------------------------------------------------------------------- 
    61  
    62       IF(lwp) THEN 
    63          WRITE(numout,*) 
    64          WRITE(numout,*) '   ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients' 
    65          WRITE(numout,*) '   ~~~~~~~' 
    66       ENDIF 
    67  
     62      ! 
     63      IF(lwp) WRITE(numout,*) 
     64      IF(lwp) WRITE(numout,*) '   ldf_c1d : set a given profile to eddy mixing coefficients' 
     65      ! 
    6866      ! initialization of the profile 
     67      zratio = 0.25_wp           ! surface/bottom ratio 
    6968      zh =  500._wp              ! depth    of the inflection point [m] 
    7069      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m] 
    7170      !                          ! associated coefficient           [-] 
    72       zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) ) 
     71      zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) ) 
    7372      ! 
    7473      ! 
     
    7675      ! 
    7776      CASE( 'DYN' )                     ! T- and F-points 
    78          DO jk = 1, jpk                      ! pah1 at T-point 
    79             pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk) 
     77         DO jk = jpkm1, 1, -1                ! pah1 at T-point 
     78            pah1(:,:,jk) = pahs1(:,:) * (  zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) )  ) 
    8079         END DO 
    81          DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.) 
     80         DO jk = jpkm1, 1, -1                ! pah2 at F-point (zdep2 is an approximation in zps-coord.) 
    8281            DO jj = 1, jpjm1 
    83                DO ji = 1, fs_jpim1 
    84                   zdep2 = (  gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk)   & 
    85                      &     + gdept_n(ji,jj  ,jk) + gdept_n(ji+1,jj  ,jk)  ) * 0.25_wp 
    86                   pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk) 
     82               DO ji = 1, jpim1 
     83                  zdep2 = (  gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk)   & 
     84                     &     + gdept_0(ji,jj  ,jk) + gdept_0(ji+1,jj  ,jk)  ) * r1_4 
     85                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) 
    8786               END DO 
    8887            END DO 
     
    9190         ! 
    9291      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) 
    93          DO jk = 1, jpk 
     92         DO jk = jpkm1, 1, -1 
    9493            DO jj = 1, jpjm1 
    95                DO ji = 1, fs_jpim1 
    96                   zdep1 = (  gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk)  ) * 0.5_wp 
    97                   zdep2 = (  gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk)  ) * 0.5_wp 
    98                   pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk) 
    99                   pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk) 
     94               DO ji = 1, jpim1 
     95                  zdep1 = (  gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk)  ) * 0.5_wp 
     96                  zdep2 = (  gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk)  ) * 0.5_wp 
     97                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) 
     98                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) 
    10099               END DO 
    101100            END DO 
     
    104103         CALL lbc_lnk_multi( pah1, 'U', 1. , pah2, 'V', 1. )    
    105104         ! 
     105      CASE DEFAULT                        ! error 
     106         CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 
    106107      END SELECT 
    107108      ! 
     
    109110 
    110111 
    111    SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 ) 
     112   SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 ) 
    112113      !!---------------------------------------------------------------------- 
    113114      !!                  ***  ROUTINE ldf_c2d  *** 
     
    124125      !!             DYN      pah1, pah2 defined at T- and F-points 
    125126      !!---------------------------------------------------------------------- 
    126       CHARACTER(len=3)                , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers 
    127       CHARACTER(len=3)                , INTENT(in   ) ::   cd_op        ! operator: LAPlacian BiLaPlacian 
    128       REAL(wp)                        , INTENT(in   ) ::   pah0         ! eddy coefficient   [m2/s or m4/s] 
    129       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pah1, pah2   ! eddy coefficient   [m2/s or m4/s] 
     127      CHARACTER(len=3)          , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers 
     128      REAL(wp)                  , INTENT(in   ) ::   pUfac        ! =1/2*Uc LAPlacian BiLaPlacian 
     129      INTEGER                   , INTENT(in   ) ::   knn          ! characteristic velocity   [m/s] 
     130      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pah1, pah2   ! eddy coefficients         [m2/s or m4/s] 
    130131      ! 
    131132      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    132       REAL(wp) ::   za00, zd_max, zemax1, zemax2   ! local scalar 
     133      INTEGER  ::   inn          ! local integer 
    133134      !!---------------------------------------------------------------------- 
    134135      ! 
    135       zd_max = ra * rad       ! = 1 degree at the equator in meters 
     136      IF(lwp) WRITE(numout,*) 
     137      IF(lwp) WRITE(numout,*) '   ldf_c2d :   aht = Ufac * max(e1,e2)   with Ufac = ', pUfac, ' m/s' 
    136138      ! 
    137       IF(lwp) THEN 
    138          WRITE(numout,*) 
    139          WRITE(numout,*) '   ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) ' 
    140          WRITE(numout,*) '   ~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)' 
    141          WRITE(numout,*) 
    142       ENDIF 
    143139      ! 
    144       SELECT CASE( cd_type )        !==  surface values  ==!  (depending on DYN/TRA) 
     140      SELECT CASE( cd_type )        !==  surface values  ==!  (chosen grid point function of DYN or TRA) 
    145141      ! 
    146       CASE( 'DYN' )                     ! T- and F-points 
    147          IF( cd_op == 'LAP' ) THEN            ! laplacian operator 
    148             IF(lwp) WRITE(numout,*) '              momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 
    149             za00 = pah0 / zd_max 
    150             DO jj = 1, jpj  
    151                DO ji = 1, jpi  
    152                   zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 
    153                   zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 
    154                   pah1(ji,jj,1) = za00 * zemax1 
    155                   pah2(ji,jj,1) = za00 * zemax2 
    156                END DO 
     142      CASE( 'DYN' )                       ! T- and F-points 
     143         DO jj = 1, jpj 
     144            DO ji = 1, jpi  
     145               pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn 
     146               pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn 
    157147            END DO 
    158          ELSEIF( cd_op == 'BLP' ) THEN     ! bilaplacian operator 
    159             IF(lwp) WRITE(numout,*) '              momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 
    160             za00 = pah0 / ( zd_max * zd_max * zd_max ) 
    161             DO jj = 1, jpj 
    162                DO ji = 1, jpi 
    163                   zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 
    164                   zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 
    165                   pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1  
    166                   pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2  
    167                END DO 
     148         END DO 
     149      CASE( 'TRA' )                       ! U- and V-points 
     150         DO jj = 1, jpj  
     151            DO ji = 1, jpi  
     152               pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 
     153               pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn 
    168154            END DO 
    169          ELSE                                ! NO diffusion/viscosity 
    170             CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) 
    171          ENDIF 
    172          !                                !  deeper values  (LAP and BLP cases) 
    173          DO jk = 2, jpk 
    174             pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk)  
    175             pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk)  
    176155         END DO 
    177          ! 
    178       CASE( 'TRA' )                     ! U- and V-points (approximation in zps-coord.) 
    179          IF( cd_op == 'LAP' ) THEN            ! laplacian operator 
    180             IF(lwp) WRITE(numout,*) '              tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 
    181             za00 = pah0 / zd_max 
    182             DO jj = 1, jpj  
    183                DO ji = 1, jpi  
    184                   zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 
    185                   zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 
    186                   pah1(ji,jj,1) = za00 * zemax1 
    187                   pah2(ji,jj,1) = za00 * zemax2 
    188                END DO 
    189             END DO 
    190          ELSEIF( cd_op == 'BLP' ) THEN      ! bilaplacian operator (NB: square root of the coeff) 
    191             IF(lwp) WRITE(numout,*) '              tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 
    192             za00 = pah0 / ( zd_max * zd_max * zd_max ) 
    193             DO jj = 1, jpj 
    194                DO ji = 1, jpi 
    195                   zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1)  
    196                   zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1)  
    197                   pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1  
    198                   pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2  
    199                END DO 
    200             END DO 
    201          ELSE                                ! NO diffusion/viscosity 
    202             CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) 
    203          ENDIF 
    204          !                                !  deeper values  (LAP and BLP cases) 
    205          DO jk = 2, jpk 
    206             pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk)  
    207             pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk)  
    208          END DO 
    209          ! 
     156      CASE DEFAULT                        ! error 
     157         CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 
    210158      END SELECT 
     159      !                             !==  deeper values = surface one  ==!  (except jpk) 
     160      DO jk = 2, jpkm1 
     161         pah1(:,:,jk) = pah1(:,:,1) 
     162         pah2(:,:,jk) = pah2(:,:,1) 
     163      END DO 
    211164      ! 
    212165   END SUBROUTINE ldf_c2d 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90

    r9169 r9490  
    1717   USE dom_oce         ! ocean space and time domain  
    1818   USE phycst          ! physical constants 
     19   USE ldfslp          ! lateral diffusion: slopes of mixing orientation 
    1920   USE ldfc1d_c2d      ! lateral diffusion: 1D and 2D cases 
    2021   ! 
     
    3132   PUBLIC   ldf_dyn        ! called by step.F90 
    3233 
    33    !                                                !!* Namelist namdyn_ldf : lateral mixing on momentum * 
     34   !                                    !!* Namelist namdyn_ldf : lateral mixing on momentum * 
    3435   LOGICAL , PUBLIC ::   ln_dynldf_NONE  !: No operator (i.e. no explicit diffusion) 
    3536   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator 
     
    3738   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction 
    3839   LOGICAL , PUBLIC ::   ln_dynldf_hor   !: horizontal (geopotential) direction 
    39    LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction 
     40!  LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction                        (see ldfslp) 
    4041   INTEGER , PUBLIC ::   nn_ahm_ijk_t    !: choice of time & space variations of the lateral eddy viscosity coef. 
    41    REAL(wp), PUBLIC ::   rn_ahm_0        !: lateral laplacian eddy viscosity            [m2/s] 
    42    REAL(wp), PUBLIC ::   rn_ahm_b        !: lateral laplacian background eddy viscosity [m2/s] 
    43    REAL(wp), PUBLIC ::   rn_bhm_0        !: lateral bilaplacian eddy viscosity          [m4/s] 
    44                                          !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 
    45                                          !! will be computed. 
    46    REAL(wp), PUBLIC ::   rn_csmc         !: Smagorinsky constant of proportionality  
    47    REAL(wp), PUBLIC ::   rn_minfac       !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 
    48    REAL(wp), PUBLIC ::   rn_maxfac       !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 
    49  
    50    LOGICAL , PUBLIC ::   l_ldfdyn_time   !: flag for time variation of the lateral eddy viscosity coef. 
    51  
    52    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy diffusivity coef. at U- and V-points   [m2/s or m4/s] 
     42   !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case)  
     43   !                                           !                             bht = 1/12 Ud*Ld^3 (blp case) 
     44   REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s] 
     45   REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m] 
     46   !                                        ! Smagorinsky viscosity  (nn_ahm_ijk_t = 32)  
     47   REAL(wp), PUBLIC ::   rn_csmc               !: Smagorinsky constant of proportionality  
     48   REAL(wp), PUBLIC ::   rn_minfac             !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 
     49   REAL(wp), PUBLIC ::   rn_maxfac             !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 
     50   !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 
     51   REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s] 
     52 
     53   !                                    !!* Parameter to control the type of lateral viscous operator 
     54   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator 
     55   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00                       !: without operator (i.e. no lateral viscous trend) 
     56   !                          !!      laplacian     !    bilaplacian    ! 
     57   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator 
     58   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator 
     59   ! 
     60   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
     61   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef. 
     62 
     63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 
    5364   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
    5465   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only) 
    5566   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2            
    5667 
    57    REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12 
     68   REAL(wp) ::   r1_2    = 0.5_wp            ! =1/2 
    5869   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4 
    5970   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8 
     71   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12 
    6072   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 ) 
    6173 
     
    92104      !!                                                           or  L^4|D|/8  bilaplacian operator ) 
    93105      !!---------------------------------------------------------------------- 
    94       INTEGER  ::   ji, jj, jk        ! dummy loop indices 
    95       INTEGER  ::   ierr, inum, ios   ! local integer 
    96       REAL(wp) ::   zah0              ! local scalar 
    97       ! 
     106      INTEGER  ::   ji, jj, jk                     ! dummy loop indices 
     107      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer 
     108      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar 
     109      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s) 
     110      !! 
    98111      NAMELIST/namdyn_ldf/ ln_dynldf_NONE, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator 
    99          &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso ,   &   ! acting direction of the operator 
    100          &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 ,   &   ! lateral eddy coefficient 
    101          &                 rn_csmc      , rn_minfac, rn_maxfac                 ! Smagorinsky settings 
     112         &                 ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator 
     113         &                 nn_ahm_ijk_t  , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient 
     114         &                 rn_csmc       , rn_minfac    , rn_maxfac            ! Smagorinsky settings 
    102115      !!---------------------------------------------------------------------- 
    103116      ! 
     
    129142         WRITE(numout,*) '      coefficients :' 
    130143         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t 
    131          WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0      = ', rn_ahm_0, ' m2/s' 
    132          WRITE(numout,*) '         background viscosity (iso case)      rn_ahm_b      = ', rn_ahm_b, ' m2/s' 
    133          WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_bhm_0      = ', rn_bhm_0, ' m4/s' 
     144         WRITE(numout,*) '         lateral viscous velocity  (if cst)      rn_Uv      = ', rn_Uv, ' m/s' 
     145         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m' 
     146         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s' 
     147         ! 
    134148         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :' 
    135149         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc 
    136          WRITE(numout,*) '         factor multiplier for theorectical lower limit for ' 
    137          WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_minfac    = ', rn_minfac 
    138          WRITE(numout,*) '         factor multiplier for theorectical lower upper for ' 
    139          WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_maxfac    = ', rn_maxfac 
     150         WRITE(numout,*) '         factor multiplier for eddy visc.' 
     151         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac 
     152         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac 
    140153      ENDIF 
    141154 
    142       !                                ! Parameter control 
     155      ! 
     156      !           !==  type of lateral operator used  ==!   (set nldf_dyn) 
     157      !           !=====================================! 
     158      ! 
     159      nldf_dyn = np_ERROR 
     160      ioptio = 0 
     161      IF( ln_dynldf_NONE ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF 
     162      IF( ln_dynldf_lap  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
     163      IF( ln_dynldf_blp  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
     164      IF( ioptio /= 1    )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
     165      ! 
     166      IF(.NOT.ln_dynldf_NONE ) THEN    !==  direction ==>> type of operator  ==! 
     167         ioptio = 0 
     168         IF( ln_dynldf_lev )   ioptio = ioptio + 1 
     169         IF( ln_dynldf_hor )   ioptio = ioptio + 1 
     170         IF( ln_dynldf_iso )   ioptio = ioptio + 1 
     171         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 
     172         ! 
     173         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals 
     174         ierr = 0 
     175         IF( ln_dynldf_lap ) THEN         ! laplacian operator 
     176            IF( ln_zco ) THEN                   ! z-coordinate 
     177               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation) 
     178               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation) 
     179               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation) 
     180            ENDIF 
     181            IF( ln_zps ) THEN                   ! z-coordinate with partial step 
     182               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level              (no rotation) 
     183               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level              (no rotation) 
     184               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation) 
     185            ENDIF 
     186            IF( ln_sco ) THEN                   ! s-coordinate 
     187               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation) 
     188               IF ( ln_dynldf_hor )   nldf_dyn = np_lap_i   ! horizontal             (   rotation) 
     189               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation) 
     190            ENDIF 
     191         ENDIF 
     192         ! 
     193         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator 
     194            IF( ln_zco ) THEN                   ! z-coordinate 
     195               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation) 
     196               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation) 
     197               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation) 
     198            ENDIF 
     199            IF( ln_zps ) THEN                   ! z-coordinate with partial step 
     200               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation) 
     201               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level              (no rotation) 
     202               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation) 
     203            ENDIF 
     204            IF( ln_sco ) THEN                   ! s-coordinate 
     205               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation) 
     206               IF( ln_dynldf_hor )   ierr = 2            ! horizontal             (   rotation) 
     207               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation) 
     208            ENDIF 
     209         ENDIF 
     210         ! 
     211         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 
     212         ! 
     213         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes 
     214         ! 
     215      ENDIF 
     216      ! 
     217      IF(lwp) THEN 
     218         WRITE(numout,*) 
     219         SELECT CASE( nldf_dyn ) 
     220         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity' 
     221         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator' 
     222         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background' 
     223         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator' 
     224         END SELECT 
     225         WRITE(numout,*) 
     226      ENDIF 
     227       
     228      ! 
     229      !           !==  Space/time variation of eddy coefficients  ==! 
     230      !           !=================================================! 
     231      ! 
     232      l_ldfdyn_time = .FALSE.                ! no time variation except in case defined below 
     233      ! 
    143234      IF( ln_dynldf_NONE ) THEN 
    144235         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated' 
    145          l_ldfdyn_time = .FALSE. 
    146236         RETURN 
    147       ENDIF 
    148       ! 
    149       IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN     ! iso-neutral bilaplacian not implemented 
    150          CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' )  
    151       ENDIF 
    152  
    153       ! ... Space/Time variation of eddy coefficients 
    154       !                                               ! allocate the ahm arrays 
    155       ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 
    156       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 
    157       ! 
    158       ahmt(:,:,jpk) = 0._wp                           ! last level always 0   
    159       ahmf(:,:,jpk) = 0._wp 
    160       ! 
    161       !                                               ! value of eddy mixing coef. 
    162       IF    ( ln_dynldf_lap ) THEN   ;   zah0 =      rn_ahm_0         !   laplacian operator 
    163       ELSEIF( ln_dynldf_blp ) THEN   ;   zah0 = ABS( rn_bhm_0 )       ! bilaplacian operator 
    164       ELSE                                                                  ! NO viscous  operator 
    165          CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 
    166       ENDIF 
    167       ! 
    168       l_ldfdyn_time = .FALSE.                         ! no time variation except in case defined below 
    169       ! 
    170       IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN     ! only if a lateral diffusion operator is used 
    171          ! 
    172          SELECT CASE(  nn_ahm_ijk_t  )                ! Specification of space time variations of ahmt, ahmf 
     237         ! 
     238      ELSE                                   !==  a lateral diffusion operator is used  ==! 
     239         ! 
     240         !                                         ! allocate the ahm arrays 
     241         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 
     242         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 
     243         ! 
     244         ahmt(:,:,jpk) = 0._wp                     ! last level always 0   
     245         ahmf(:,:,jpk) = 0._wp 
     246         ! 
     247         !                                         ! value of lap/blp eddy mixing coef. 
     248         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian 
     249         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian 
     250         ENDIF 
     251         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient 
     252         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator) 
     253         ! 
     254         SELECT CASE(  nn_ahm_ijk_t  )             !* Specification of space-time variations of ahmt, ahmf 
    173255         ! 
    174256         CASE(   0  )      !==  constant  ==! 
    175             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = constant ' 
    176             ahmt(:,:,:) = zah0 * tmask(:,:,:) 
    177             ahmf(:,:,:) = zah0 * fmask(:,:,:) 
     257            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity. = constant = ', zah0, cl_Units 
     258            ahmt(:,:,1:jpkm1) = zah0 
     259            ahmf(:,:,1:jpkm1) = zah0 
    178260            ! 
    179261         CASE(  10  )      !==  fixed profile  ==! 
    180             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F( depth )' 
    181             ahmt(:,:,1) = zah0 * tmask(:,:,1)                      ! constant surface value 
    182             ahmf(:,:,1) = zah0 * fmask(:,:,1) 
    183             CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 
     262            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( depth )' 
     263            IF(lwp) WRITE(numout,*) '           surface viscous coef. = constant = ', zah0, cl_Units 
     264            ahmt(:,:,1) = zah0                        ! constant surface value 
     265            ahmf(:,:,1) = zah0 
     266            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 
    184267            ! 
    185268         CASE ( -20 )      !== fixed horizontal shape read in file  ==! 
    186             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file' 
     269            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file' 
    187270            CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 
    188271            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 
    189272            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 
    190273            CALL iom_close( inum ) 
    191 !!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ??? 
    192 !!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 
    193 !!              better:  check that the max is <=1  i.e. it is a shape from 0 to 1, not a coef that has physical dimension 
    194274            DO jk = 2, jpkm1 
    195                ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk) 
    196                ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk) 
     275               ahmt(:,:,jk) = ahmt(:,:,1) 
     276               ahmf(:,:,jk) = ahmf(:,:,1) 
    197277            END DO 
    198278            ! 
    199279         CASE(  20  )      !== fixed horizontal shape  ==! 
    200             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 
    201             IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor 
    202             IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor^3 
     280            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 
     281            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Lv = Max(e1,e2)' 
     282            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)' 
     283            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn 
    203284            ! 
    204285         CASE( -30  )      !== fixed 3D shape read in file  ==! 
    205             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 
     286            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file' 
    206287            CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 
    207288            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 
    208289            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 
    209290            CALL iom_close( inum ) 
    210 !!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ???? 
    211 !!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 
    212             DO jk = 1, jpkm1 
    213                ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk) 
    214                ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk) 
    215             END DO 
    216291            ! 
    217292         CASE(  30  )       !==  fixed 3D shape  ==! 
    218             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F( latitude, longitude, depth )' 
    219             IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor 
    220             IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor 
    221             !                                                    ! reduction with depth 
    222             CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 
     293            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth )' 
     294            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Ld = Max(e1,e2)' 
     295            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)' 
     296            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn 
     297            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )  ! reduction with depth 
    223298            ! 
    224299         CASE(  31  )       !==  time varying 3D field  ==! 
    225             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F( latitude, longitude, depth , time )' 
    226             IF(lwp) WRITE(numout,*) '              proportional to the velocity : |u|e/12 or |u|e^3/12' 
     300            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )' 
     301            IF(lwp) WRITE(numout,*) '           proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)' 
    227302            ! 
    228303            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90 
    229304            ! 
    230305         CASE(  32  )       !==  time varying 3D field  ==! 
    231             IF(lwp) WRITE(numout,*) '   ==>>>   momentum mixing coef. = F( latitude, longitude, depth , time )' 
    232             IF(lwp) WRITE(numout,*) '              proportional to the local deformation rate and gridscale (Smagorinsky)' 
    233             IF(lwp) WRITE(numout,*) '                                                                : L^2|D| or L^4|D|/8' 
     306            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )' 
     307            IF(lwp) WRITE(numout,*) '           proportional to the local deformation rate and gridscale (Smagorinsky)' 
    234308            ! 
    235309            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90 
    236310            ! 
    237             ! allocate arrays used in ldf_dyn.  
     311            !                          ! allocate arrays used in ldf_dyn.  
    238312            ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) ,  esqf(jpi,jpj) , STAT=ierr ) 
    239313            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 
    240314            ! 
    241             ! Set local gridscale values 
    242             DO jj = 2, jpjm1 
     315            DO jj = 2, jpjm1           ! Set local gridscale values 
    243316               DO ji = fs_2, fs_jpim1 
    244317                  esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2  
     
    251324         END SELECT 
    252325         ! 
    253          IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN       ! bilapcian and no time variation: 
    254             ahmt(:,:,:) = SQRT( ahmt(:,:,:) )                     ! take the square root of the coefficient 
    255             ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 
     326         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation  
     327            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only) 
     328               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
     329               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
     330            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask) 
     331               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 
     332               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 
     333            ENDIF 
    256334         ENDIF 
    257335         ! 
     
    341419                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  & 
    342420                          &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk) 
    343                      dtensq(ji,jj) = zdb*zdb 
     421                     dtensq(ji,jj) = zdb * zdb 
    344422                  END DO 
    345423               END DO 
     
    351429                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  & 
    352430                          &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk) 
    353                      dshesq(ji,jj) = zdb*zdb 
     431                     dshesq(ji,jj) = zdb * zdb 
    354432                  END DO 
    355433               END DO 
     
    385463         ! 
    386464         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 
    387                                        !                      = sqrt( A_lap_smag L^2/8 ) 
    388                                        ! stability limits already applied to laplacian values 
    389                                        ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 
    390             ! 
     465            !                          !                      = sqrt( A_lap_smag L^2/8 ) 
     466            !                          ! stability limits already applied to laplacian values 
     467            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4 
    391468            DO jk = 1, jpkm1 
    392                ! 
    393469               DO jj = 2, jpjm1 
    394470                  DO ji = fs_2, fs_jpim1 
    395                      ! 
    396                      ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 
    397                      ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 
    398                      ! 
     471                     ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 
     472                     ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 
    399473                  END DO 
    400474               END DO 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r9124 r9490  
    2222   USE oce            ! ocean dynamics and tracers 
    2323   USE dom_oce        ! ocean space and time domain 
    24    USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
     24!   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    2525   USE phycst         ! physical constants 
    2626   USE zdfmxl         ! mixed layer depth 
     
    4343   LOGICAL , PUBLIC ::   l_ldfslp = .FALSE.     !: slopes flag 
    4444 
    45    LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.       !: iso-neutral direction 
    46    LOGICAL , PUBLIC ::   ln_traldf_triad = .FALSE.      !: griffies triad scheme 
    47  
    48    LOGICAL , PUBLIC ::   ln_triad_iso    = .FALSE.      !: pure horizontal mixing in ML 
    49    LOGICAL , PUBLIC ::   ln_botmix_triad = .FALSE.      !: mixing on bottom 
    50    REAL(wp), PUBLIC ::   rn_sw_triad     = 1._wp        !: =1 switching triads ; =0 all four triads used  
    51    REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp      !: slope limit 
     45   LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.       !: iso-neutral direction                           (nam_traldf namelist) 
     46   LOGICAL , PUBLIC ::   ln_traldf_triad = .FALSE.      !: griffies triad scheme                           (nam_traldf namelist) 
     47   LOGICAL , PUBLIC ::   ln_dynldf_iso                  !: iso-neutral direction                           (nam_dynldf namelist) 
     48 
     49   LOGICAL , PUBLIC ::   ln_triad_iso    = .FALSE.      !: pure horizontal mixing in ML                    (nam_traldf namelist) 
     50   LOGICAL , PUBLIC ::   ln_botmix_triad = .FALSE.      !: mixing on bottom                                (nam_traldf namelist) 
     51   REAL(wp), PUBLIC ::   rn_sw_triad     = 1._wp        !: =1 switching triads ; =0 all four triads used   (nam_traldf namelist) 
     52   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp      !: slope limit                                     (nam_traldf namelist) 
    5253 
    5354   LOGICAL , PUBLIC ::   l_grad_zps = .FALSE.           !: special treatment for Horz Tgradients w partial steps (triad operator) 
     
    749750      ! 
    750751      IF( ln_traldf_triad ) THEN        ! Griffies operator : triad of slopes 
    751          IF(lwp) WRITE(numout,*) '              Griffies (triad) operator initialisation' 
     752         IF(lwp) WRITE(numout,*) '   ==>>>   triad) operator (Griffies)' 
    752753         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) ,     & 
    753754            &      triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1) ,     & 
     
    757758         ! 
    758759      ELSE                             ! Madec operator : slopes at u-, v-, and w-points 
    759          IF(lwp) WRITE(numout,*) '              Madec operator initialisation' 
     760         IF(lwp) WRITE(numout,*) '   ==>>>   iso operator (Madec)' 
    760761         ALLOCATE( omlmask(jpi,jpj,jpk) ,                                                                        & 
    761762            &      uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) ,     & 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90

    r9169 r9490  
    5050   LOGICAL , PUBLIC ::   ln_traldf_hor       !: horizontal (geopotential) direction 
    5151!  LOGICAL , PUBLIC ::   ln_traldf_iso       !: iso-neutral direction                    (see ldfslp) 
     52   !                                    != iso-neutral options =! 
    5253!  LOGICAL , PUBLIC ::   ln_traldf_triad     !: griffies triad scheme                    (see ldfslp) 
    5354   LOGICAL , PUBLIC ::   ln_traldf_msc       !: Method of Stabilizing Correction  
     
    5859   !                                    !=  Coefficients =! 
    5960   INTEGER , PUBLIC ::   nn_aht_ijk_t        !: choice of time & space variations of the lateral eddy diffusivity coef. 
    60    REAL(wp), PUBLIC ::   rn_aht_0            !:   laplacian lateral eddy diffusivity [m2/s] 
    61    REAL(wp), PUBLIC ::   rn_bht_0            !: bilaplacian lateral eddy diffusivity [m4/s] 
    62  
    63    !                                   !!* Namelist namtra_ldfeiv : eddy induced velocity param. * 
     61   !                                            !  time invariant coefficients:  aht_0 = 1/2  Ud*Ld   (lap case)  
     62   !                                            !                                bht_0 = 1/12 Ud*Ld^3 (blp case) 
     63   REAL(wp), PUBLIC ::      rn_Ud               !: lateral diffusive velocity  [m/s] 
     64   REAL(wp), PUBLIC ::      rn_Ld               !: lateral diffusive length    [m] 
     65 
     66   !                                   !!* Namelist namtra_eiv : eddy induced velocity param. * 
    6467   !                                    != Use/diagnose eiv =! 
    6568   LOGICAL , PUBLIC ::   ln_ldfeiv           !: eddy induced velocity flag 
     
    6770   !                                    != Coefficients =! 
    6871   INTEGER , PUBLIC ::   nn_aei_ijk_t        !: choice of time/space variation of the eiv coeff. 
    69    REAL(wp), PUBLIC ::   rn_aeiv_0           !: eddy induced velocity coefficient [m2/s] 
     72   REAL(wp), PUBLIC ::      rn_Ue               !: lateral diffusive velocity  [m/s] 
     73   REAL(wp), PUBLIC ::      rn_Le               !: lateral diffusive length    [m] 
    7074    
     75   !                                  ! Flag to control the type of lateral diffusive operator 
     76   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in specification of lateral diffusion 
     77   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend) 
     78   !                          !!      laplacian     !    bilaplacian    ! 
     79   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator 
     80   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator 
     81   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator 
     82 
     83   INTEGER , PUBLIC ::   nldf_tra      = 0         !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 
    7184   LOGICAL , PUBLIC ::   l_ldftra_time = .FALSE.   !: flag for time variation of the lateral eddy diffusivity coef. 
    72    LOGICAL , PUBLIC ::   l_ldfeiv_time = .FALSE.   ! flag for time variation of the eiv coef. 
     85   LOGICAL , PUBLIC ::   l_ldfeiv_time = .FALSE.   !: flag for time variation of the eiv coef. 
    7386 
    7487   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahtu, ahtv   !: eddy diffusivity coef. at U- and V-points   [m2/s] 
    7588   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aeiu, aeiv   !: eddy induced velocity coeff.                [m2/s] 
    7689 
     90   REAL(wp) ::   aht0, aei0               ! constant eddy coefficients (deduced from namelist values)                     [m2/s] 
     91   REAL(wp) ::   r1_2  = 0.5_wp           ! =1/2 
    7792   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4 
    7893   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12 
     
    108123      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file 
    109124      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10) 
    110       !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator 
    111       !!                                                          or |u|e^3/12 bilaplacian operator ) 
     125      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  1/2  |u|e     laplacian operator 
     126      !!                                                           or 1/12 |u|e^3 bilaplacian operator ) 
    112127      !!              * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init  
    113128      !!             
    114       !! ** action  : ahtu, ahtv initialized once for all or l_ldftra_time set to true 
    115       !!              aeiu, aeiv initialized once for all or l_ldfeiv_time set to true 
    116       !!---------------------------------------------------------------------- 
    117       INTEGER  ::   jk                ! dummy loop indices 
    118       INTEGER  ::   ierr, inum, ios   ! local integer 
    119       REAL(wp) ::   zah0              ! local scalar 
     129      !! ** action  : ahtu, ahtv initialized one for all or l_ldftra_time set to true 
     130      !!              aeiu, aeiv initialized one for all or l_ldfeiv_time set to true 
     131      !!---------------------------------------------------------------------- 
     132      INTEGER  ::   jk                             ! dummy loop indices 
     133      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer 
     134      REAL(wp) ::   zah_max, zUfac                 !   -      - 
     135      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s) 
    120136      !! 
    121137      NAMELIST/namtra_ldf/ ln_traldf_NONE, ln_traldf_lap  , ln_traldf_blp  ,  &   ! type of operator 
     
    123139         &                 ln_traldf_iso , ln_traldf_msc  ,  rn_slpmax     ,  &   ! option for iso-neutral operator 
    124140         &                 ln_triad_iso  , ln_botmix_triad, rn_sw_triad    ,  &   ! option for triad operator 
    125          &                 rn_aht_0      , rn_bht_0       , nn_aht_ijk_t          ! lateral eddy coefficient 
     141         &                 nn_aht_ijk_t  , rn_Ud          , rn_Ld                 ! lateral eddy coefficient 
    126142      !!---------------------------------------------------------------------- 
    127143      ! 
    128144      IF(lwp) THEN                      ! control print 
    129145         WRITE(numout,*) 
    130          WRITE(numout,*) 'ldf_tra_init : lateral tracer physics' 
     146         WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion' 
    131147         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    132148      ENDIF 
     149       
    133150      ! 
    134151      !  Choice of lateral tracer physics 
     
    154171         WRITE(numout,*) '         iso-neutral Madec operator              ln_traldf_iso   = ', ln_traldf_iso 
    155172         WRITE(numout,*) '         iso-neutral triad operator              ln_traldf_triad = ', ln_traldf_triad 
    156          WRITE(numout,*) '            iso-neutral (Method of Stab. Corr.)  ln_traldf_msc   = ', ln_traldf_msc 
     173         WRITE(numout,*) '            use the Method of Stab. Correction   ln_traldf_msc   = ', ln_traldf_msc 
    157174         WRITE(numout,*) '            maximum isoppycnal slope             rn_slpmax       = ', rn_slpmax 
    158175         WRITE(numout,*) '            pure lateral mixing in ML            ln_triad_iso    = ', ln_triad_iso 
     
    160177         WRITE(numout,*) '            lateral mixing on bottom             ln_botmix_triad = ', ln_botmix_triad 
    161178         WRITE(numout,*) '      coefficients :' 
    162          WRITE(numout,*) '         lateral eddy diffusivity   (lap case)   rn_aht_0        = ', rn_aht_0 
    163          WRITE(numout,*) '         lateral eddy diffusivity (bilap case)   rn_bht_0        = ', rn_bht_0 
    164179         WRITE(numout,*) '         type of time-space variation            nn_aht_ijk_t    = ', nn_aht_ijk_t 
    165       ENDIF 
    166       ! 
    167       !                                ! Parameter control 
    168       ! 
    169       IF( ln_traldf_NONE ) THEN 
    170          IF(lwp) WRITE(numout,*) '   ==>>>   No diffusive operator selected. ahtu and ahtv are not allocated' 
    171          l_ldftra_time = .FALSE. 
    172          RETURN 
    173       ENDIF 
     180         WRITE(numout,*) '            lateral diffusive velocity (if cst)  rn_Ud           = ', rn_Ud, ' m/s' 
     181         WRITE(numout,*) '            lateral diffusive length   (if cst)  rn_Ld           = ', rn_Ld, ' m' 
     182      ENDIF 
     183      ! 
     184      ! 
     185      ! Operator and its acting direction   (set nldf_tra)   
     186      ! ================================= 
     187      ! 
     188      nldf_tra = np_ERROR 
     189      ioptio   = 0 
     190      IF( ln_traldf_NONE ) THEN   ;   nldf_tra = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF 
     191      IF( ln_traldf_lap  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
     192      IF( ln_traldf_blp  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
     193      IF( ioptio /=  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
     194      ! 
     195      IF( .NOT.ln_traldf_NONE ) THEN   !==  direction ==>> type of operator  ==! 
     196         ioptio = 0 
     197         IF( ln_traldf_lev )   ioptio = ioptio + 1 
     198         IF( ln_traldf_hor )   ioptio = ioptio + 1 
     199         IF( ln_traldf_iso )   ioptio = ioptio + 1 
     200         IF( ioptio /=  1  )   CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 
     201         ! 
     202         !                                ! defined the type of lateral diffusion from ln_traldf_... logicals 
     203         ierr = 0 
     204         IF ( ln_traldf_lap ) THEN        ! laplacian operator 
     205            IF ( ln_zco ) THEN                  ! z-coordinate 
     206               IF ( ln_traldf_lev   )   nldf_tra = np_lap     ! iso-level = horizontal (no rotation) 
     207               IF ( ln_traldf_hor   )   nldf_tra = np_lap     ! iso-level = horizontal (no rotation) 
     208               IF ( ln_traldf_iso   )   nldf_tra = np_lap_i   ! iso-neutral: standard  (   rotation) 
     209               IF ( ln_traldf_triad )   nldf_tra = np_lap_it  ! iso-neutral: triad     (   rotation) 
     210            ENDIF 
     211            IF ( ln_zps ) THEN                  ! z-coordinate with partial step 
     212               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed  
     213               IF ( ln_traldf_hor   )   nldf_tra = np_lap     ! horizontal             (no rotation) 
     214               IF ( ln_traldf_iso   )   nldf_tra = np_lap_i   ! iso-neutral: standard     (rotation) 
     215               IF ( ln_traldf_triad )   nldf_tra = np_lap_it  ! iso-neutral: triad        (rotation) 
     216            ENDIF 
     217            IF ( ln_sco ) THEN                  ! s-coordinate 
     218               IF ( ln_traldf_lev   )   nldf_tra = np_lap     ! iso-level              (no rotation) 
     219               IF ( ln_traldf_hor   )   nldf_tra = np_lap_i   ! horizontal             (   rotation) 
     220               IF ( ln_traldf_iso   )   nldf_tra = np_lap_i   ! iso-neutral: standard  (   rotation) 
     221               IF ( ln_traldf_triad )   nldf_tra = np_lap_it  ! iso-neutral: triad     (   rotation) 
     222            ENDIF 
     223         ENDIF 
     224         ! 
     225         IF( ln_traldf_blp ) THEN         ! bilaplacian operator 
     226            IF ( ln_zco ) THEN                  ! z-coordinate 
     227               IF ( ln_traldf_lev   )   nldf_tra = np_blp     ! iso-level = horizontal (no rotation) 
     228               IF ( ln_traldf_hor   )   nldf_tra = np_blp     ! iso-level = horizontal (no rotation) 
     229               IF ( ln_traldf_iso   )   nldf_tra = np_blp_i   ! iso-neutral: standard  (   rotation) 
     230               IF ( ln_traldf_triad )   nldf_tra = np_blp_it  ! iso-neutral: triad     (   rotation) 
     231            ENDIF 
     232            IF ( ln_zps ) THEN                  ! z-coordinate with partial step 
     233               IF ( ln_traldf_lev   )   ierr     = 1          ! iso-level not allowed  
     234               IF ( ln_traldf_hor   )   nldf_tra = np_blp     ! horizontal             (no rotation) 
     235               IF ( ln_traldf_iso   )   nldf_tra = np_blp_i   ! iso-neutral: standard  (   rotation) 
     236               IF ( ln_traldf_triad )   nldf_tra = np_blp_it  ! iso-neutral: triad     (   rotation) 
     237            ENDIF 
     238            IF ( ln_sco ) THEN                  ! s-coordinate 
     239               IF ( ln_traldf_lev   )   nldf_tra = np_blp     ! iso-level              (no rotation) 
     240               IF ( ln_traldf_hor   )   nldf_tra = np_blp_it  ! horizontal             (   rotation) 
     241               IF ( ln_traldf_iso   )   nldf_tra = np_blp_i   ! iso-neutral: standard  (   rotation) 
     242               IF ( ln_traldf_triad )   nldf_tra = np_blp_it  ! iso-neutral: triad     (   rotation) 
     243            ENDIF 
     244         ENDIF 
     245         IF ( ierr == 1 )   CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 
     246      ENDIF 
     247      ! 
     248      IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                & 
     249           &            CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 
     250      IF( ln_isfcav .AND. ln_traldf_triad ) & 
     251           &            CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 
     252           ! 
     253      IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 
     254         & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
    174255      ! 
    175256      IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN     ! iso-neutral bilaplacian need MSC 
     
    177258      ENDIF 
    178259      ! 
     260      IF(lwp) THEN 
     261         WRITE(numout,*) 
     262         SELECT CASE( nldf_tra ) 
     263         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral diffusion' 
     264         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   laplacian iso-level operator' 
     265         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   Rotated laplacian operator (standard)' 
     266         CASE( np_lap_it )   ;   WRITE(numout,*) '   ==>>>   Rotated laplacian operator (triad)' 
     267         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   bilaplacian iso-level operator' 
     268         CASE( np_blp_i  )   ;   WRITE(numout,*) '   ==>>>   Rotated bilaplacian operator (standard)' 
     269         CASE( np_blp_it )   ;   WRITE(numout,*) '   ==>>>   Rotated bilaplacian operator (triad)' 
     270         END SELECT 
     271         WRITE(numout,*) 
     272      ENDIF 
     273 
     274      ! 
    179275      !  Space/time variation of eddy coefficients  
    180276      ! =========================================== 
    181       !                                               ! allocate the aht arrays 
    182       ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 
    183       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 
    184       ! 
    185       ahtu(:,:,jpk) = 0._wp                           ! last level always 0   
    186       ahtv(:,:,jpk) = 0._wp 
    187       ! 
    188       !                                               ! value of eddy mixing coef. 
    189       IF    ( ln_traldf_lap ) THEN   ;   zah0 =      rn_aht_0        !   laplacian operator 
    190       ELSEIF( ln_traldf_blp ) THEN   ;   zah0 = ABS( rn_bht_0 )      ! bilaplacian operator 
    191       ENDIF 
    192       ! 
    193       l_ldftra_time = .FALSE.                         ! no time variation except in case defined below 
    194       ! 
    195       IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN     ! only if a lateral diffusion operator is used 
    196          ! 
    197          SELECT CASE(  nn_aht_ijk_t  )                   ! Specification of space time variations of ehtu, ahtv 
     277      ! 
     278      l_ldftra_time = .FALSE.                ! no time variation except in case defined below 
     279      ! 
     280      IF( ln_traldf_NONE ) THEN              !== no explicit diffusive operator  ==! 
     281         ! 
     282         IF(lwp) WRITE(numout,*) '   ==>>>   No diffusive operator selected. ahtu and ahtv are not allocated' 
     283         RETURN 
     284         ! 
     285      ELSE                                   !==  a lateral diffusion operator is used  ==! 
     286         ! 
     287         !                                         ! allocate the aht arrays 
     288         ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 
     289         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 
     290         ! 
     291         ahtu(:,:,jpk) = 0._wp                     ! last level always 0   
     292         ahtv(:,:,jpk) = 0._wp 
     293         !. 
     294         !                                         ! value of lap/blp eddy mixing coef. 
     295         IF(     ln_traldf_lap ) THEN   ;   zUfac = r1_2 *rn_Ud   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian 
     296         ELSEIF( ln_traldf_blp ) THEN   ;   zUfac = r1_12*rn_Ud   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian 
     297         ENDIF 
     298         aht0    = zUfac *    rn_Ld**inn              ! mixing coefficient 
     299         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator for e1=1 degree) 
     300         ! 
     301         ! 
     302         SELECT CASE(  nn_aht_ijk_t  )             !* Specification of space-time variations of ahtu, ahtv 
    198303         ! 
    199304         CASE(   0  )      !==  constant  ==! 
    200             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = constant = ', rn_aht_0 
    201             ahtu(:,:,:) = zah0 * umask(:,:,:) 
    202             ahtv(:,:,:) = zah0 * vmask(:,:,:) 
     305            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = constant = ', aht0, cl_Units 
     306            ahtu(:,:,1:jpkm1) = aht0 
     307            ahtv(:,:,1:jpkm1) = aht0 
    203308            ! 
    204309         CASE(  10  )      !==  fixed profile  ==! 
    205             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( depth )' 
    206             ahtu(:,:,1) = zah0 * umask(:,:,1)                      ! constant surface value 
    207             ahtv(:,:,1) = zah0 * vmask(:,:,1) 
    208             CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 
    209             ! 
    210          CASE ( -20 )      !== fixed horizontal shape read in file  ==! 
    211             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file' 
     310            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F( depth )' 
     311            IF(lwp) WRITE(numout,*) '           surface eddy diffusivity = constant = ', aht0, cl_Units 
     312            ahtu(:,:,1) = aht0                        ! constant surface value 
     313            ahtv(:,:,1) = aht0 
     314            CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 
     315            ! 
     316         CASE ( -20 )      !== fixed horizontal shape and magnitude read in file  ==! 
     317            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file' 
    212318            CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 
    213319            CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) ) 
     
    215321            CALL iom_close( inum ) 
    216322            DO jk = 2, jpkm1 
    217                ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 
    218                ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) 
     323               ahtu(:,:,jk) = ahtu(:,:,1) 
     324               ahtv(:,:,jk) = ahtv(:,:,1) 
    219325            END DO 
    220326            ! 
    221327         CASE(  20  )      !== fixed horizontal shape  ==! 
    222             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 
    223             IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor 
    224             IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor 
     328            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 
     329            IF(lwp) WRITE(numout,*) '           using a fixed diffusive velocity = ', rn_Ud,' m/s   and   Ld = Max(e1,e2)' 
     330            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)' 
     331            CALL ldf_c2d( 'TRA', zUfac      , inn        , ahtu, ahtv )    ! value proportional to scale factor^inn 
    225332            ! 
    226333         CASE(  21  )      !==  time varying 2D field  ==! 
    227             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, time )' 
    228             IF(lwp) WRITE(numout,*) '                               = F( growth rate of baroclinic instability )' 
    229             IF(lwp) WRITE(numout,*) '                               min value = 0.1 * rn_aht_0' 
    230             IF(lwp) WRITE(numout,*) '                               max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)' 
    231             IF(lwp) WRITE(numout,*) '                               increased to rn_aht_0 within 20N-20S' 
     334            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F( latitude, longitude, time )' 
     335            IF(lwp) WRITE(numout,*) '                            = F( growth rate of baroclinic instability )' 
     336            IF(lwp) WRITE(numout,*) '            min value = 0.2 * aht0 (with aht0= 1/2 rn_Ud*rn_Ld)' 
     337            IF(lwp) WRITE(numout,*) '            max value =       aei0 (with aei0=1/2 rn_Ue*Le  increased to aht0 within 20N-20S' 
    232338            ! 
    233339            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90 
    234340            ! 
    235             IF( ln_traldf_blp ) THEN 
    236                CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' ) 
    237             ENDIF 
     341            IF( ln_traldf_blp )   CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)',   & 
     342               &                                 '              incompatible with bilaplacian operator' ) 
    238343            ! 
    239344         CASE( -30  )      !== fixed 3D shape read in file  ==! 
    240             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j,k) read in eddy_diffusivity.nc file' 
     345            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file' 
    241346            CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 
    242347            CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu ) 
    243348            CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv ) 
    244349            CALL iom_close( inum ) 
    245             DO jk = 1, jpkm1 
    246                ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk) 
    247                ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 
    248             END DO 
    249350            ! 
    250351         CASE(  30  )      !==  fixed 3D shape  ==! 
    251             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth )' 
    252             IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor 
    253             IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor 
    254             !                                                    ! reduction with depth 
    255             CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 
     352            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F( latitude, longitude, depth )' 
     353            IF(lwp) WRITE(numout,*) '           using a fixed diffusive velocity = ', rn_Ud,' m/s   and   Ld = Max(e1,e2)' 
     354            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)' 
     355            CALL ldf_c2d( 'TRA', zUfac      , inn        , ahtu, ahtv )    ! surface value proportional to scale factor^inn 
     356            CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )    ! reduction with depth 
    256357            ! 
    257358         CASE(  31  )      !==  time varying 3D field  ==! 
    258             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth , time )' 
    259             IF(lwp) WRITE(numout,*) '                                 proportional to the velocity : |u|e/12 or |u|e^3/12' 
     359            IF(lwp) WRITE(numout,*) '   ==>>>   eddy diffusivity = F( latitude, longitude, depth , time )' 
     360            IF(lwp) WRITE(numout,*) '           proportional to the velocity : 1/2 |u|e or 1/12 |u|e^3' 
    260361            ! 
    261362            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90 
     
    265366         END SELECT 
    266367         ! 
    267          IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN 
    268             ahtu(:,:,:) = SQRT( ahtu(:,:,:) ) 
    269             ahtv(:,:,:) = SQRT( ahtv(:,:,:) ) 
     368         IF( .NOT.l_ldftra_time ) THEN             !* No time variation  
     369            IF(     ln_traldf_lap ) THEN                 !   laplacian operator (mask only) 
     370               ahtu(:,:,1:jpkm1) =       ahtu(:,:,1:jpkm1)   * umask(:,:,1:jpkm1) 
     371               ahtv(:,:,1:jpkm1) =       ahtv(:,:,1:jpkm1)   * vmask(:,:,1:jpkm1) 
     372            ELSEIF( ln_traldf_blp ) THEN                 ! bilaplacian operator (square root + mask) 
     373               ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1) 
     374               ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1) 
     375            ENDIF 
    270376         ENDIF 
    271377         ! 
     
    281387      !! ** Purpose :   update at kt the tracer lateral mixing coeff. (aht and aeiv) 
    282388      !! 
    283       !! ** Method  :   time varying eddy diffusivity coefficients: 
     389      !! ** Method  : * time varying eddy diffusivity coefficients: 
    284390      !! 
    285391      !!    nn_aei_ijk_t = 21    aeiu, aeiv = F(i,j,  t) = F(growth rate of baroclinic instability) 
     
    290396      !!                                                                     or |u|e^3/12 bilaplacian operator ) 
    291397      !! 
     398      !!              * time varying EIV coefficients: call to ldf_eiv routine 
     399      !! 
    292400      !! ** action  :   ahtu, ahtv   update at each time step    
    293401      !!                aeiu, aeiv      -       -     -    -   (if ln_ldfeiv=T)  
     
    296404      ! 
    297405      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    298       REAL(wp) ::   zaht, zahf, zaht_min, z1_f20       ! local scalar 
     406      REAL(wp) ::   zaht, zahf, zaht_min, zDaht, z1_f20   ! local scalar 
    299407      !!---------------------------------------------------------------------- 
    300408      ! 
    301409      IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN       ! eddy induced velocity coefficients 
    302410         !                                ! =F(growth rate of baroclinic instability) 
    303          !                                ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S 
    304          CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv ) 
     411         !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S 
     412         CALL ldf_eiv( kt, aei0, aeiu, aeiv ) 
    305413      ENDIF 
    306414      ! 
     
    308416      ! 
    309417      CASE(  21  )       !==  time varying 2D field  ==!   = F( growth rate of baroclinic instability ) 
    310          !                                             !   min value rn_aht_0 / 10  
    311          !                                             !   max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21) 
    312          !                                             !   increase to rn_aht_0 within 20N-20S 
     418         !                                             !   min value 0.2*aht0 
     419         !                                             !   max value aht0 (aei0 if nn_aei_ijk_t=21) 
     420         !                                             !   increase to aht0 within 20N-20S 
    313421         IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN   ! use the already computed aei. 
    314422            ahtu(:,:,1) = aeiu(:,:,1) 
    315423            ahtv(:,:,1) = aeiv(:,:,1) 
    316424         ELSE                                            ! compute aht.  
    317             CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv ) 
     425            CALL ldf_eiv( kt, aht0, ahtu, ahtv ) 
    318426         ENDIF 
    319427         ! 
    320          z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )      ! 1 / ff(20 degrees)    
    321          zaht_min = 0.2_wp * rn_aht_0                                      ! minimum value for aht 
     428         z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )   ! 1 / ff(20 degrees)    
     429         zaht_min = 0.2_wp * aht0                                       ! minimum value for aht 
     430         zDaht    = aht0 - zaht_min                                       
    322431         DO jj = 1, jpj 
    323432            DO ji = 1, jpi 
    324433               !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 
    325434               !!     ==>>>   The Coriolis value is identical for t- & u_points, and for v- and f-points 
    326                zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 
    327                zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 
    328                ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  ) * umask(ji,jj,1)     ! min value zaht_min 
    329                ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  ) * vmask(ji,jj,1)     ! increase within 20S-20N 
    330             END DO 
    331          END DO 
    332          DO jk = 2, jpkm1                             ! deeper value = surface value 
     435               zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 
     436               zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 
     437               ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  )     ! min value zaht_min 
     438               ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  )     ! increase within 20S-20N 
     439            END DO 
     440         END DO 
     441         DO jk = 1, jpkm1                             ! deeper value = surface value + mask for all levels 
    333442            ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 
    334443            ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) 
     
    338447         IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12 
    339448            DO jk = 1, jpkm1 
    340                ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 
     449               ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12   ! n.b. ub,vb are masked 
    341450               ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 
    342451            END DO 
     
    355464      CALL iom_put( "ahtv_3d", ahtv(:,:,:) )   ! 3D      v-eddy diffusivity coeff. 
    356465      ! 
    357 !!gm  : THE IF below is to be checked (comes from Seb) 
    358466      IF( ln_ldfeiv ) THEN 
    359467        CALL iom_put( "aeiu_2d", aeiu(:,:,1) )   ! surface u-EIV coeff. 
     
    372480      !! ** Purpose :   initialization of the eiv coeff. from namelist choices. 
    373481      !! 
    374       !! ** Method : 
    375       !! 
    376       !! ** Action :   aeiu , aeiv   : EIV coeff. at u- & v-points 
    377       !!               l_ldfeiv_time : =T if EIV coefficients vary with time 
    378       !!---------------------------------------------------------------------- 
    379       INTEGER  ::   jk                ! dummy loop indices 
    380       INTEGER  ::   ierr, inum, ios   ! local integer 
    381       ! 
    382       NAMELIST/namtra_ldfeiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &    ! eddy induced velocity (eiv) 
    383          &                    nn_aei_ijk_t, rn_aeiv_0             ! eiv  coefficient 
     482      !! ** Method  :   the eiv diffusivity coef. specification depends on: 
     483      !!    nn_aei_ijk_t  =  0 => = constant 
     484      !!                  ! 
     485      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth  
     486      !!                  ! 
     487      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file 
     488      !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 
     489      !!                  = 21    = F(i,j,t) = F(growth rate of baroclinic instability) 
     490      !!                  ! 
     491      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file 
     492      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10) 
     493      !! 
     494      !! ** Action  :   aeiu , aeiv   :  initialized one for all or l_ldftra_time set to true 
     495      !!                l_ldfeiv_time : =T if EIV coefficients vary with time 
     496      !!---------------------------------------------------------------------- 
     497      INTEGER  ::   jk                     ! dummy loop indices 
     498      INTEGER  ::   ierr, inum, ios, inn   ! local integer 
     499      REAL(wp) ::   zah_max, zUfac         !   -   scalar 
     500      !! 
     501      NAMELIST/namtra_eiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &   ! eddy induced velocity (eiv) 
     502         &                 nn_aei_ijk_t, rn_Ue, rn_Le         ! eiv  coefficient 
    384503      !!---------------------------------------------------------------------- 
    385504      ! 
     
    390509      ENDIF 
    391510      ! 
    392       REWIND( numnam_ref )              ! Namelist namtra_ldfeiv in reference namelist : eddy induced velocity param. 
    393       READ  ( numnam_ref, namtra_ldfeiv, IOSTAT = ios, ERR = 901) 
    394 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_ldfeiv in reference namelist', lwp ) 
    395       ! 
    396       REWIND( numnam_cfg )              ! Namelist namtra_ldfeiv in configuration namelist : eddy induced velocity param. 
    397       READ  ( numnam_cfg, namtra_ldfeiv, IOSTAT = ios, ERR = 902 ) 
    398 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_ldfeiv in configuration namelist', lwp ) 
    399       IF(lwm)  WRITE ( numond, namtra_ldfeiv ) 
     511      REWIND( numnam_ref )              ! Namelist namtra_eiv in reference namelist : eddy induced velocity param. 
     512      READ  ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 
     513901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp ) 
     514      ! 
     515      REWIND( numnam_cfg )              ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param. 
     516      READ  ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 
     517902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist', lwp ) 
     518      IF(lwm)  WRITE ( numond, namtra_eiv ) 
    400519 
    401520      IF(lwp) THEN                      ! control print 
    402          WRITE(numout,*) '   Namelist namtra_ldfeiv : ' 
    403          WRITE(numout,*) '      Eddy Induced Velocity (eiv) param.      ln_ldfeiv     = ', ln_ldfeiv 
    404          WRITE(numout,*) '      eiv streamfunction & velocity diag.     ln_ldfeiv_dia = ', ln_ldfeiv_dia 
    405          WRITE(numout,*) '      eddy induced velocity coef.             rn_aeiv_0     = ', rn_aeiv_0 
    406          WRITE(numout,*) '      type of time-space variation            nn_aei_ijk_t  = ', nn_aei_ijk_t 
     521         WRITE(numout,*) '   Namelist namtra_eiv : ' 
     522         WRITE(numout,*) '      Eddy Induced Velocity (eiv) param.         ln_ldfeiv     = ', ln_ldfeiv 
     523         WRITE(numout,*) '      eiv streamfunction & velocity diag.        ln_ldfeiv_dia = ', ln_ldfeiv_dia 
     524         WRITE(numout,*) '      coefficients :' 
     525         WRITE(numout,*) '         type of time-space variation            nn_aei_ijk_t  = ', nn_aht_ijk_t 
     526         WRITE(numout,*) '         lateral diffusive velocity (if cst)     rn_Ue         = ', rn_Ue, ' m/s' 
     527         WRITE(numout,*) '         lateral diffusive length   (if cst)     rn_Le         = ', rn_Le, ' m' 
    407528         WRITE(numout,*) 
    408529      ENDIF 
    409530      ! 
    410       IF( ln_ldfeiv .AND. ln_traldf_blp )   CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 
    411  
    412       !                                 ! Parameter control 
    413       l_ldfeiv_time = .FALSE.     
    414       ! 
    415       IF( ln_ldfeiv ) THEN                         ! allocate the aei arrays 
     531      l_ldfeiv_time = .FALSE.       ! no time variation except in case defined below 
     532      ! 
     533      ! 
     534      IF( .NOT.ln_ldfeiv ) THEN     !== Parametrization not used  ==! 
     535         ! 
     536         IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity param is NOT used' 
     537         ln_ldfeiv_dia = .FALSE. 
     538         ! 
     539      ELSE                          !== use the parametrization  ==! 
     540         ! 
     541         IF(lwp) WRITE(numout,*) '   ==>>>   use eddy induced velocity parametrization' 
     542         IF(lwp) WRITE(numout,*) 
     543         ! 
     544         IF( ln_traldf_blp )   CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 
     545         ! 
     546         !                                != allocate the aei arrays 
    416547         ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 
    417548         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays') 
    418549         ! 
    419          SELECT CASE( nn_aei_ijk_t )               ! Specification of space time variations of eaiu, aeiv 
    420          ! 
    421          CASE(   0  )      !==  constant  ==! 
    422             IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = constant = ', rn_aeiv_0 
    423             aeiu(:,:,:) = rn_aeiv_0 
    424             aeiv(:,:,:) = rn_aeiv_0 
    425             ! 
    426          CASE(  10  )      !==  fixed profile  ==! 
     550         !                                != Specification of space-time variations of eaiu, aeiv 
     551         ! 
     552         aeiu(:,:,jpk) = 0._wp               ! last level always 0   
     553         aeiv(:,:,jpk) = 0._wp 
     554         !                                   ! value of EIV coef. (laplacian operator) 
     555         zUfac = r1_2 *rn_Ue                    ! velocity factor 
     556         inn = 1                                ! L-exponent 
     557         aei0    = zUfac *    rn_Le**inn        ! mixing coefficient 
     558         zah_max = zUfac * (ra*rad)**inn        ! maximum reachable coefficient (value at the Equator) 
     559 
     560         SELECT CASE( nn_aei_ijk_t )         !* Specification of space-time variations 
     561         ! 
     562         CASE(   0  )                        !--  constant  --! 
     563            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = constant = ', aei0, ' m2/s' 
     564            aeiu(:,:,1:jpkm1) = aei0 
     565            aeiv(:,:,1:jpkm1) = aei0 
     566            ! 
     567         CASE(  10  )                        !--  fixed profile  --! 
    427568            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( depth )' 
    428             aeiu(:,:,1) = rn_aeiv_0                                ! constant surface value 
    429             aeiv(:,:,1) = rn_aeiv_0 
    430             CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 
    431             ! 
    432          CASE ( -20 )      !== fixed horizontal shape read in file  ==! 
    433             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 
     569            IF(lwp) WRITE(numout,*) '           surface eddy diffusivity = constant = ', aht0, ' m2/s' 
     570            aeiu(:,:,1) = aei0                  ! constant surface value 
     571            aeiv(:,:,1) = aei0 
     572            CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 
     573            ! 
     574         CASE ( -20 )                        !--  fixed horizontal shape read in file  --! 
     575            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 
    434576            CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 
    435577            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) ) 
    436578            CALL iom_get  ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) ) 
    437579            CALL iom_close( inum ) 
    438             DO jk = 2, jpk 
     580            DO jk = 2, jpkm1 
    439581               aeiu(:,:,jk) = aeiu(:,:,1) 
    440582               aeiv(:,:,jk) = aeiv(:,:,1) 
    441583            END DO 
    442584            ! 
    443          CASE(  20  )      !== fixed horizontal shape  ==! 
    444             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)' 
    445             CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor 
    446             ! 
    447          CASE(  21  )       !==  time varying 2D field  ==! 
    448             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, time )' 
    449             IF(lwp) WRITE(numout,*) '                               = F( growth rate of baroclinic instability )' 
     585         CASE(  20  )                        !--  fixed horizontal shape  --! 
     586            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( e1, e2 )' 
     587            IF(lwp) WRITE(numout,*) '           using a fixed diffusive velocity = ', rn_Ue, ' m/s   and   Le = Max(e1,e2)' 
     588            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s   for e1=1°)' 
     589            CALL ldf_c2d( 'TRA', zUfac      , inn        , aeiu, aeiv )    ! value proportional to scale factor^inn 
     590            ! 
     591         CASE(  21  )                        !--  time varying 2D field  --! 
     592            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( latitude, longitude, time )' 
     593            IF(lwp) WRITE(numout,*) '                                       = F( growth rate of baroclinic instability )' 
     594            IF(lwp) WRITE(numout,*) '           maximum allowed value: aei0 = ', aei0, ' m2/s' 
    450595            ! 
    451596            l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90 
    452597            ! 
    453          CASE( -30  )      !== fixed 3D shape read in file  ==! 
    454             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 
     598         CASE( -30  )                        !-- fixed 3D shape read in file  --! 
     599            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 
    455600            CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 
    456601            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu ) 
     
    458603            CALL iom_close( inum ) 
    459604            ! 
    460          CASE(  30  )       !==  fixed 3D shape  ==! 
    461             IF(lwp) WRITE(numout,*) '   ==>>>   tracer mixing coef. = F( latitude, longitude, depth )' 
    462             CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor 
    463             !                                                 ! reduction with depth 
    464             CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 
     605         CASE(  30  )                        !--  fixed 3D shape  --! 
     606            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( latitude, longitude, depth )' 
     607            CALL ldf_c2d( 'TRA', zUfac      , inn        , aeiu, aeiv )    ! surface value proportional to scale factor^inn 
     608            CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )    ! reduction with depth 
    465609            ! 
    466610         CASE DEFAULT 
     
    468612         END SELECT 
    469613         ! 
    470       ELSE 
    471           IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity param is NOT used neither diagnosed' 
    472           ln_ldfeiv_dia = .FALSE. 
     614         IF( .NOT.l_ldfeiv_time ) THEN             !* mask if No time variation  
     615            DO jk = 1, jpkm1 
     616               aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) 
     617               ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 
     618            END DO 
     619         ENDIF 
     620         ! 
    473621      ENDIF 
    474622      !                     
     
    493641      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    494642      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei    ! local scalars 
    495       REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zross, zaeiw   ! 2D workspace 
    496       !!---------------------------------------------------------------------- 
    497       ! 
    498       zn   (:,:) = 0._wp      ! Local initialization 
    499       zhw  (:,:) = 5._wp 
    500       zah  (:,:) = 0._wp 
    501       zross(:,:) = 0._wp 
     643      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zaeiw   ! 2D workspace 
     644      !!---------------------------------------------------------------------- 
     645      ! 
     646      zn (:,:) = 0._wp        ! Local initialization 
     647      zhw(:,:) = 5._wp 
     648      zah(:,:) = 0._wp 
     649      zRo(:,:) = 0._wp 
    502650      !                       ! Compute lateral diffusive coefficient at T-point 
    503651      IF( ln_traldf_triad ) THEN 
     
    538686            END DO 
    539687         END DO 
    540       END IF 
     688      ENDIF 
    541689 
    542690      DO jj = 2, jpjm1 
    543691         DO ji = fs_2, fs_jpim1   ! vector opt. 
    544692            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 
    545             ! Rossby radius at w-point taken < 40km and  > 2km 
    546             zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 ) 
     693            ! Rossby radius at w-point taken betwenn 2 km and  40km 
     694            zRo(ji,jj) = MAX(  2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 
    547695            ! Compute aeiw by multiplying Ro^2 and T^-1 
    548             zaeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
     696            zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
    549697         END DO 
    550698      END DO 
     
    554702      DO jj = 2, jpjm1 
    555703         DO ji = fs_2, fs_jpim1   ! vector opt. 
    556             zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)       ! tropical decrease 
     704            zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease 
    557705            zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0 
    558706         END DO 
     
    620768         DO jj = 1, jpjm1 
    621769            DO ji = 1, fs_jpim1   ! vector opt. 
    622                zpsi_uw(ji,jj,jk) = - 0.25_wp * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   & 
    623                   &                                       * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * umask(ji,jj,jk) 
    624                zpsi_vw(ji,jj,jk) = - 0.25_wp * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   & 
    625                   &                                       * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * vmask(ji,jj,jk) 
     770               zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   & 
     771                  &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * umask(ji,jj,jk) 
     772               zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   & 
     773                  &                                    * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * vmask(ji,jj,jk) 
    626774            END DO 
    627775         END DO 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r9023 r9490  
    1010   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    1111   !!---------------------------------------------------------------------- 
    12  
    13    !! * Modules used 
    14    USE par_kind, ONLY : &         ! Precision variables 
    15       & wp 
    16    USE in_out_manager             ! I/O manager 
    17    USE obs_inter_sup              ! Interpolation support 
    18    USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt 
    19       & obs_int_h2d, & 
    20       & obs_int_h2d_init 
    21    USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint 
    22       & obs_avg_h2d, & 
    23       & obs_avg_h2d_init, & 
    24       & obs_max_fpsize 
    25    USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    26       & obs_int_z1d,    & 
    27       & obs_int_z1d_spl 
    28    USE obs_const,  ONLY :    &    ! Obs fill value 
    29       & obfillflt 
    30    USE dom_oce,       ONLY : & 
    31       & glamt, glamf, & 
    32       & gphit, gphif 
    33    USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    34       & ctl_warn, ctl_stop 
    35    USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
    36       & sbc_dcy, nday_qsr 
    37    USE obs_grid,      ONLY : &  
    38       & obs_level_search      
     12   USE obs_inter_sup                                        ! Interpolation support 
     13   USE obs_inter_h2d, ONLY : obs_int_h2d, obs_int_h2d_init  ! Horizontal interpolation to the obs pt 
     14   USE obs_averg_h2d, ONLY : obs_avg_h2d, obs_avg_h2d_init, obs_max_fpsize    ! Horizontal averaging to the obs footprint 
     15   USE obs_inter_z1d, ONLY : obs_int_z1d, obs_int_z1d_spl   ! Vertical interpolation to the obs pt 
     16   USE obs_const    , ONLY : obfillflt                      ! Obs fill value 
     17   USE dom_oce,       ONLY :   glamt, glamf, gphit, gphif   ! lat/lon of ocean grid-points 
     18   USE lib_mpp,       ONLY :   ctl_warn, ctl_stop           ! Warning and stopping routines 
     19   USE sbcdcy,        ONLY :   sbc_dcy, nday_qsr            ! For calculation of where it is night-time 
     20   USE obs_grid,      ONLY :   obs_level_search      
     21   ! 
     22   USE par_kind     , ONLY :   wp   ! Precision variables 
     23   USE in_out_manager               ! I/O manager 
    3924 
    4025   IMPLICIT NONE 
    41  
    42    !! * Routine accessibility 
    4326   PRIVATE 
    4427 
    45    PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
    46       &   obs_surf_opt     ! Compute the model counterpart of surface obs 
    47  
    48    INTEGER, PARAMETER, PUBLIC :: & 
    49       & imaxavtypes = 20   ! Max number of daily avgd obs types 
     28   PUBLIC   obs_prof_opt   !: Compute the model counterpart of profile obs 
     29   PUBLIC   obs_surf_opt   !: Compute the model counterpart of surface obs 
     30 
     31   INTEGER, PARAMETER, PUBLIC ::   imaxavtypes = 20   !: Max number of daily avgd obs types 
    5032 
    5133   !!---------------------------------------------------------------------- 
    52    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     34   !! NEMO/OPA 4.0 , NEMO Consortium (2018) 
    5335   !! $Id$ 
    5436   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5537   !!---------------------------------------------------------------------- 
    56  
    5738CONTAINS 
    58  
    5939 
    6040   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
     
    6444      &                     plam1, plam2, pphi1, pphi2,           & 
    6545      &                     k1dint, k2dint, kdailyavtypes ) 
    66  
    6746      !!----------------------------------------------------------------------- 
    68       !! 
    6947      !!                     ***  ROUTINE obs_pro_opt  *** 
    7048      !! 
     
    11492      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    11593      !!----------------------------------------------------------------------- 
    116  
    117       !! * Modules used 
    11894      USE obs_profiles_def ! Definition of storage space for profile obs. 
    11995 
    12096      IMPLICIT NONE 
    12197 
    122       !! * Arguments 
    123       TYPE(obs_prof), INTENT(INOUT) :: & 
    124          & prodatqc                  ! Subset of profile data passing QC 
    125       INTEGER, INTENT(IN) :: kt      ! Time step 
    126       INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
    127       INTEGER, INTENT(IN) :: kpj 
    128       INTEGER, INTENT(IN) :: kpk 
    129       INTEGER, INTENT(IN) :: kit000  ! Number of the first time step 
    130                                      !   (kit000-1 = restart time) 
    131       INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header) 
    132       INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header) 
    133       INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 
    134       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    135          & pvar1,    &               ! Model field 1 
    136          & pvar2,    &               ! Model field 2 
    137          & pmask1,   &               ! Land-sea mask 1 
    138          & pmask2                    ! Land-sea mask 2 
    139       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    140          & plam1,    &               ! Model longitudes for variable 1 
    141          & plam2,    &               ! Model longitudes for variable 2 
    142          & pphi1,    &               ! Model latitudes for variable 1 
    143          & pphi2                     ! Model latitudes for variable 2 
    144       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    145          & pgdept, &                 ! Model array of depth T levels  
    146          & pgdepw                    ! Model array of depth W levels  
    147       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    148          & kdailyavtypes             ! Types for daily averages 
     98      TYPE(obs_prof), INTENT(inout) ::   prodatqc        ! Subset of profile data passing QC 
     99      INTEGER       , INTENT(in   ) ::   kt              ! Time step 
     100      INTEGER       , INTENT(in   ) ::   kpi, kpj, kpk   ! Model grid parameters 
     101      INTEGER       , INTENT(in   ) ::   kit000          ! Number of the first time step (kit000-1 = restart time) 
     102      INTEGER       , INTENT(in   ) ::   k1dint          ! Vertical interpolation type (see header) 
     103      INTEGER       , INTENT(in   ) ::   k2dint          ! Horizontal interpolation type (see header) 
     104      INTEGER       , INTENT(in   ) ::   kdaystp         ! Number of time steps per day 
     105      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pvar1 , pvar2    ! Model field     1 and 2 
     106      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pmask1, pmask2   ! Land-sea mask   1 and 2 
     107      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   plam1 , plam2    ! Model longitude 1 and 2 
     108      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   pphi1 , pphi2    ! Model latitudes 1 and 2 
     109      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pgdept, pgdepw   ! depth of T and W levels  
     110      INTEGER, DIMENSION(imaxavtypes), OPTIONAL ::   kdailyavtypes             ! Types for daily averages 
    149111 
    150112      !! * Local declarations 
     
    706668      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    707669      !!----------------------------------------------------------------------- 
    708  
    709       !! * Modules used 
    710670      USE obs_surf_def  ! Definition of storage space for surface observations 
    711671 
    712672      IMPLICIT NONE 
    713673 
    714       !! * Arguments 
    715674      TYPE(obs_surf), INTENT(INOUT) :: & 
    716675         & surfdataqc                  ! Subset of surface data passing QC 
     
    866825         DO ji = 0, imaxifp 
    867826            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
    868              
     827            ! 
    869828            !Deal with wrap around in longitude 
    870829            IF ( imodi < 1      ) imodi = imodi + jpiglo 
    871830            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
    872              
     831            ! 
    873832            DO jj = 0, imaxjfp 
    874833               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     
    877836               IF ( imodj < 1      ) imodj = 1 
    878837               IF ( imodj > jpjglo ) imodj = jpjglo 
    879  
     838               ! 
    880839               igrdip1(ji+1,jj+1,iobs) = imodi 
    881840               igrdjp1(ji+1,jj+1,iobs) = imodj 
    882                 
     841               ! 
    883842               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
    884843                  igrdi(ji,jj,iobs) = imodi 
    885844                  igrdj(ji,jj,iobs) = imodj 
    886845               ENDIF 
    887                 
     846               ! 
    888847            END DO 
    889848         END DO 
     
    1010969            & ) 
    1011970      ENDIF 
    1012  
     971      ! 
    1013972      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
    1014  
     973      ! 
    1015974   END SUBROUTINE obs_surf_opt 
    1016975 
     976   !!====================================================================== 
    1017977END MODULE obs_oper 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r9023 r9490  
    108108      imin0 = ( nn_time0 - ihou0 * 100 ) 
    109109 
    110       icycle = no     ! Assimilation cycle 
     110      icycle = nn_no     ! Assimilation cycle 
    111111 
    112112      ! Diagnotics counters for various failures. 
     
    339339      imin0 = ( nn_time0 - ihou0 * 100 ) 
    340340 
    341       icycle = no     ! Assimilation cycle 
     341      icycle = nn_no     ! Assimilation cycle 
    342342 
    343343      ! Diagnotics counters for various failures. 
    344344 
    345       iotdobs  = 0 
    346       igrdobs  = 0 
     345      iotdobs   = 0 
     346      igrdobs   = 0 
    347347      iosdv1obs = 0 
    348348      iosdv2obs = 0 
     
    884884      !!        !  2007-01  (K. Mogensen)  Original 
    885885      !!---------------------------------------------------------------------- 
    886       !! * Arguments 
    887886      INTEGER, INTENT(IN) :: kobsno        ! Number of observations 
    888887      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: & 
     
    924923      !! ** Action  :  
    925924      !!    
    926       !! History : 
    927       !!        !  2007-03  (A. Weaver, K. Mogensen)  Original 
    928       !!        !  2007-06  (K. Mogensen et al) Reject obs. near land. 
    929       !!---------------------------------------------------------------------- 
    930       !! * Modules used 
    931  
    932       !! * Arguments 
    933       INTEGER, INTENT(IN) :: kobsno    ! Total number of observations 
    934       INTEGER, INTENT(IN) :: kpi       ! Number of grid points in (i,j) 
    935       INTEGER, INTENT(IN) :: kpj 
    936       INTEGER, DIMENSION(kobsno), INTENT(IN) :: & 
    937          & kobsi, &           ! Observation (i,j) coordinates 
    938          & kobsj 
    939       REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: & 
    940          & pobslam, &         ! Observation (lon,lat) coordinates 
    941          & pobsphi 
    942       REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: & 
    943          & plam, pphi         ! Model (lon,lat) coordinates 
    944       REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: & 
    945          & pmask              ! Land mask array 
    946       INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    947          & kobsqc             ! Observation quality control 
    948       INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
    949       INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
    950       INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
    951       INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
    952       LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
    953       LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
    954       INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
    955  
    956       !! * Local declarations 
    957       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    958          & zgmsk              ! Grid mask 
    959  
    960       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    961          & zbmsk              ! Boundary mask 
    962       REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    963       REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    964          & zglam, &           ! Model longitude at grid points 
    965          & zgphi              ! Model latitude at grid points 
    966       INTEGER, DIMENSION(2,2,kobsno) :: & 
    967          & igrdi, &           ! Grid i,j 
    968          & igrdj 
    969       LOGICAL :: lgridobs           ! Is observation on a model grid point. 
    970       INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    971       INTEGER :: jobs, ji, jj 
     925      !! History :  2007-03  (A. Weaver, K. Mogensen)  Original 
     926      !!         !  2007-06  (K. Mogensen et al) Reject obs. near land. 
     927      !!---------------------------------------------------------------------- 
     928      INTEGER , INTENT(in   )                     ::   kobsno            ! Total number of observations 
     929      INTEGER , INTENT(in   )                     ::   kpi    , kpj      ! Number of grid points in (i,j) 
     930      INTEGER , INTENT(in   ), DIMENSION(kobsno)  ::   kobsi  , kobsj    ! Observation (i,j) coordinates 
     931      REAL(wp), INTENT(in   ), DIMENSION(kobsno)  ::   pobslam, pobsphi  ! Observation (lon,lat) coordinates 
     932      REAL(wp), INTENT(in   ), DIMENSION(kpi,kpj) ::   plam   , pphi     ! Model (lon,lat) coordinates 
     933      REAL(wp), INTENT(in   ), DIMENSION(kpi,kpj) ::   pmask             ! Land mask array 
     934      INTEGER , INTENT(inout), DIMENSION(kobsno)  ::   kobsqc            ! Observation quality control 
     935      INTEGER , INTENT(inout)                     ::   kosdobs           ! Observations outside space domain 
     936      INTEGER , INTENT(inout)                     ::   klanobs           ! Observations within a model land cell 
     937      INTEGER , INTENT(inout)                     ::   knlaobs           ! Observations near land 
     938      INTEGER , INTENT(inout)                     ::   kbdyobs           ! Observations near boundary 
     939      LOGICAL , INTENT(in   )                     ::   ld_nea            ! Flag observations near land 
     940      LOGICAL , INTENT(in   )                     ::   ld_bound_reject   ! Flag observations near open boundary  
     941      INTEGER , INTENT(in   )                     ::   kqc_cutoff        ! Cutoff QC value 
     942      ! 
     943      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zgmsk          ! Grid mask 
     944      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zbmsk          ! Boundary mask 
     945      REAL(KIND=wp), DIMENSION(jpi,jpj)    ::   zbdymask 
     946      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zglam, zgphi   ! Model Lon/lat at grid points 
     947      INTEGER      , DIMENSION(2,2,kobsno) ::   igrdi, igrdj   ! Grid i,j 
     948      LOGICAL ::   lgridobs           ! Is observation on a model grid point. 
     949      INTEGER ::   iig, ijg           ! i,j of observation on model grid point. 
     950      INTEGER ::   jobs, ji, jj 
     951      !!---------------------------------------------------------------------- 
    972952       
    973953      ! Get grid point indices 
     
    11001080            ENDIF 
    11011081         ENDIF 
    1102              
     1082         ! 
    11031083      END DO 
    1104  
     1084      ! 
    11051085   END SUBROUTINE obs_coo_spc_2d 
     1086 
    11061087 
    11071088   SUBROUTINE obs_coo_spc_3d( kprofno, kobsno,  kpstart, kpend, & 
     
    11981179      INTEGER :: iig, ijg           ! i,j of observation on model grid point. 
    11991180      INTEGER :: jobs, jobsp, jk, ji, jj 
     1181      !!---------------------------------------------------------------------- 
    12001182 
    12011183      ! Get grid point indices 
     
    13591341               ENDIF 
    13601342            ENDIF 
    1361              
     1343            ! 
    13621344         END DO 
    13631345      END DO 
    1364  
     1346      ! 
    13651347   END SUBROUTINE obs_coo_spc_3d 
     1348 
    13661349 
    13671350   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 
     
    13771360      !! References : 
    13781361      !!    
    1379       !! History : 
    1380       !!        !  2007-10  (K. Mogensen) Original code 
    1381       !!---------------------------------------------------------------------- 
    1382       !! * Modules used 
    1383       !! * Arguments 
    1384       TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data 
    1385       INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value 
    1386  
    1387       !! * Local declarations 
     1362      !! History :   2007-10  (K. Mogensen) Original code 
     1363      !!---------------------------------------------------------------------- 
     1364      TYPE(obs_prof), INTENT(inout) ::   profdata     ! Profile data 
     1365      INTEGER       , INTENT(in   ) ::   kqc_cutoff   ! QC cutoff value 
     1366      ! 
    13881367      INTEGER :: jprof 
    13891368      INTEGER :: jvar 
    13901369      INTEGER :: jobs 
     1370      !!---------------------------------------------------------------------- 
    13911371       
    13921372      ! Loop over profiles 
     
    14111391 
    14121392      END DO 
    1413  
     1393      ! 
    14141394   END SUBROUTINE obs_pro_rej 
     1395 
    14151396 
    14161397   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
     
    14261407      !! References : 
    14271408      !!    
    1428       !! History : 
    1429       !!        !  2009-2  (K. Mogensen) Original code 
    1430       !!---------------------------------------------------------------------- 
    1431       !! * Modules used 
    1432       !! * Arguments 
     1409      !! History :   2009-2  (K. Mogensen) Original code 
     1410      !!---------------------------------------------------------------------- 
    14331411      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data 
    14341412      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected 
    14351413      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected 
    14361414      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
    1437  
    1438       !! * Local declarations 
     1415      ! 
    14391416      INTEGER :: jprof 
    14401417      INTEGER :: jvar 
    14411418      INTEGER :: jobs 
    1442        
    1443       ! Loop over profiles 
    1444  
    1445       DO jprof = 1, profdata%nprof 
    1446  
     1419      !!---------------------------------------------------------------------- 
     1420 
     1421      DO jprof = 1, profdata%nprof      !==  Loop over profiles  ==! 
     1422         ! 
    14471423         IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. & 
    14481424            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN 
    1449  
     1425            ! 
    14501426            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej') 
    14511427            RETURN 
    1452  
    1453          ENDIF 
    1454  
     1428            ! 
     1429         ENDIF 
     1430         ! 
    14551431         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 
    1456              
     1432             
    14571433            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
    14581434               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     
    14651441               knumu = knumu + 1 
    14661442            ENDIF 
    1467              
     1443            ! 
    14681444         END DO 
    1469              
     1445         ! 
    14701446      END DO 
    1471  
     1447      ! 
    14721448   END SUBROUTINE obs_uv_rej 
    14731449 
     1450   !!===================================================================== 
    14741451END MODULE obs_prep 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r9168 r9490  
    197197      ENDIF 
    198198      ! 
    199       IF( ln_tradmp) THEN 
     199      IF( ln_tradmp ) THEN 
    200200         !                          ! Allocate arrays 
    201201         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 
     
    212212         !    so can damp to something other than intitial conditions files? 
    213213         !!gm: In principle yes. Nevertheless, we can't anticipate demands that have never been formulated. 
    214          IF( .NOT.ln_tsd_tradmp ) THEN 
     214         IF( .NOT.ln_tsd_dmp ) THEN 
    215215            IF(lwp) WRITE(numout,*) 
    216             IF(lwp) WRITE(numout, *)  '   read T-S data not initialized, we force ln_tsd_tradmp=T' 
     216            IF(lwp) WRITE(numout, *)  '   read T-S data not initialized, we force ln_tsd_dmp=T' 
    217217            CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data 
    218218         ENDIF 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r9190 r9490  
    3737   PUBLIC   tra_ldf        ! called by step.F90  
    3838   PUBLIC   tra_ldf_init   ! called by nemogcm.F90  
    39    ! 
    40    INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 
    4139    
    4240   !! * Substitutions 
    4341#  include "vectopt_loop_substitute.h90" 
    4442   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     43   !! NEMO/OPA 4.0 , NEMO Consortium (2018) 
    4644   !! $Id$  
    4745   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    6866      ENDIF 
    6967      ! 
    70       SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
     68      SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend 
    7169      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    7270         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
     
    7674         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
    7775      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
    78          CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf ) 
     76         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra ) 
    7977      END SELECT 
    8078      ! 
     
    10199      !! ** Purpose :   Choice of the operator for the lateral tracer diffusion 
    102100      !! 
    103       !! ** Method  :   set nldf from the namtra_ldf logicals 
     101      !! ** Method  :   set nldf_tra from the namtra_ldf logicals 
    104102      !!---------------------------------------------------------------------- 
    105103      INTEGER ::   ioptio, ierr   ! temporary integers  
     
    112110         WRITE(numout,*) '   Namelist namtra_ldf: already read in ldftra module' 
    113111         WRITE(numout,*) '      see ldf_tra_init report for lateral mixing parameters' 
    114       ENDIF 
    115       !                                !==  use of lateral operator or not  ==! 
    116       nldf   = np_ERROR 
    117       ioptio = 0 
    118       IF( ln_traldf_NONE ) THEN   ;   nldf = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF 
    119       IF( ln_traldf_lap  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF 
    120       IF( ln_traldf_blp  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF 
    121       IF( ioptio /=  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
    122       ! 
    123       IF( .NOT.ln_traldf_NONE ) THEN   !==  direction ==>> type of operator  ==! 
    124          ioptio = 0 
    125          IF( ln_traldf_lev )   ioptio = ioptio + 1 
    126          IF( ln_traldf_hor )   ioptio = ioptio + 1 
    127          IF( ln_traldf_iso )   ioptio = ioptio + 1 
    128          IF( ioptio /=  1  )   CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 
     112         WRITE(numout,*) 
    129113         ! 
    130          !                                ! defined the type of lateral diffusion from ln_traldf_... logicals 
    131          ierr = 0 
    132          IF( ln_traldf_lap ) THEN         ! laplacian operator 
    133             IF ( ln_zco ) THEN               ! z-coordinate 
    134                IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
    135                IF ( ln_traldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
    136                IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
    137                IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
    138             ENDIF 
    139             IF ( ln_zps ) THEN               ! z-coordinate with partial step 
    140                IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed  
    141                IF ( ln_traldf_hor   )   nldf = np_lap     ! horizontal             (no rotation) 
    142                IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard     (rotation) 
    143                IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad        (rotation) 
    144             ENDIF 
    145             IF ( ln_sco ) THEN               ! s-coordinate 
    146                IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level              (no rotation) 
    147                IF ( ln_traldf_hor   )   nldf = np_lap_i   ! horizontal             (   rotation) 
    148                IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
    149                IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
    150             ENDIF 
    151          ENDIF 
    152          ! 
    153          IF( ln_traldf_blp ) THEN         ! bilaplacian operator 
    154             IF ( ln_zco ) THEN               ! z-coordinate 
    155                IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
    156                IF ( ln_traldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
    157                IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
    158                IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    159             ENDIF 
    160             IF ( ln_zps ) THEN               ! z-coordinate with partial step 
    161                IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed  
    162                IF ( ln_traldf_hor   )   nldf = np_blp     ! horizontal             (no rotation) 
    163                IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
    164                IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    165             ENDIF 
    166             IF ( ln_sco ) THEN               ! s-coordinate 
    167                IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level              (no rotation) 
    168                IF ( ln_traldf_hor   )   nldf = np_blp_it  ! horizontal             (   rotation) 
    169                IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
    170                IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
    171             ENDIF 
    172          ENDIF 
    173          IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 
    174       ENDIF 
    175       ! 
    176       IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                                    & 
    177            &            CALL ctl_stop( 'eddy induced velocity on tracers requires iso-neutral laplacian diffusion' ) 
    178       IF( ln_isfcav .AND. ln_traldf_triad ) & 
    179            &            CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 
    180            ! 
    181       IF(  nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 
    182          & nldf == np_blp_i .OR. nldf == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
    183       ! 
    184       IF(lwp) THEN 
    185          WRITE(numout,*) 
    186          SELECT CASE( nldf ) 
     114         SELECT CASE( nldf_tra )             ! print the choice of operator 
    187115         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral diffusion' 
    188116         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   laplacian iso-level operator' 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90

    r9124 r9490  
    3333   PUBLIC   tra_ldf_blp   ! called by traldf.F90 
    3434 
    35    !                      ! Flag to control the type of lateral diffusive operator 
    36    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in specification of lateral diffusion 
    37    INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend) 
    38    !                          !!      laplacian     !    bilaplacian    ! 
    39    INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator 
    40    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator 
    41    INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator 
    42  
    4335   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
    4436   LOGICAL  ::   l_hst   ! flag to compute heat transport 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r9436 r9490  
    234234      INTEGER  ::   ji                 ! dummy loop indices 
    235235      INTEGER  ::   ios, ilocal_comm   ! local integers 
    236       CHARACTER(len=120), DIMENSION(30) ::   cltxt, cltxt2, clnam 
     236      CHARACTER(len=120), DIMENSION(60) ::   cltxt, cltxt2, clnam 
    237237      !! 
    238238      NAMELIST/namctl/ ln_ctl   , nn_print, nn_ictls, nn_ictle,   & 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step.F90

    r9485 r9490  
    122122 
    123123      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    124       ! Ocean physics update                (ua, va, tsa used as workspace) 
     124      ! Ocean physics update 
    125125      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    126126      !  THERMODYNAMICS 
     
    208208       
    209209      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    210       ! diagnostics and outputs             (ua, va, tsa used as workspace) 
    211       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    212       IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats 
    213       IF( ln_diacfl  )   CALL dia_cfl( kstp )         ! Courant number diagnostics 
    214       IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth) 
    215       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports 
    216                          CALL dia_ar5( kstp )         ! ar5 diag 
     210      ! diagnostics and outputs 
     211      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     212      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats 
     213      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
     214      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
     215      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
     216                         CALL dia_ar5 ( kstp )        ! ar5 diag 
    217217      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    218                          CALL dia_wri( kstp )         ! ocean model: outputs 
     218                         CALL dia_wri ( kstp )        ! ocean model: outputs 
    219219      ! 
    220220      IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output 
Note: See TracChangeset for help on using the changeset viewer.