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

Changeset 14054


Ignore:
Timestamp:
2020-12-03T14:55:50+01:00 (4 years ago)
Author:
ayoung
Message:

Updated to trunk at 14052. No conflicts since last sette test at 14034. Ticket #2567.

Location:
NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC
Files:
21 deleted
109 edited
9 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r13558 r14054  
    299299!----------------------------------------------------------------------- 
    300300   ln_dynvor_een = .true.  !  energy & enstrophy scheme 
    301       nn_een_e3f = 0          ! =0   e3f = mean masked e3t divided by 4 
    302301/ 
    303302!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r13558 r14054  
    300300!----------------------------------------------------------------------- 
    301301   ln_dynvor_een = .true.  !  energy & enstrophy scheme 
    302       nn_een_e3f = 0          ! =0   e3f = mean masked e3t divided by 4 
    303302/ 
    304303!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/AMM12/EXPREF/namelist_cfg

    r13558 r14054  
    291291!----------------------------------------------------------------------- 
    292292   ln_dynvor_een = .true.  !  energy & enstrophy scheme 
    293       nn_een_e3f = 1             !  e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     293   nn_e3f_typ = 1          !  e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    294294/ 
    295295!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r14049 r14054  
    334334!----------------------------------------------------------------------- 
    335335   ln_dynvor_een = .true.  !  energy & enstrophy scheme 
    336       nn_een_e3f = 0          ! =0   e3f = mean masked e3t divided by 4 
    337336/ 
    338337!----------------------------------------------------------------------- 
    … …  
    389388      !                       !                 = 3 as =2 with distinct dissipative an mixing length scale 
    390389      nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to NIWs 
    391                                !        = 0 none ; = 1 add a tke source below the ML 
    392                                !        = 2 add a tke source just at the base of the ML 
    393                                !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T) 
     390      !                       !        = 0 none ; = 1 add a tke source below the ML 
     391      !                       !        = 2 add a tke source just at the base of the ML 
     392      !                       !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T) 
    394393      ln_mxhsw    = .false.   !  surface mixing length scale = F(wave height) 
    395394/ 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg

    r12845 r14054  
    2020! 
    2121   ln_trcdta     =  .true.  !  Initialisation from data input file (T) or not (F) 
    22    ln_trcbc      =  .false. !  Enables Boundary conditions 
     22   ln_trcbc      =  .true. !  Enables Boundary conditions 
    2323!                !           !                                           !             !         ! 
    24 !                !    name   !           title of the field              !   units     ! init    ! sbc    ! cbc    !  obc  !  
    25    sn_tracer(1)   = 'DIC     ' , 'Dissolved inorganic Concentration      ',  'mol-C/L' , .true.  , .false., .true. , .false.  
    26    sn_tracer(2)   = 'Alkalini' , 'Total Alkalinity Concentration         ',  'eq/L '   , .true.  , .false., .true. , .false. 
    27    sn_tracer(3)   = 'O2      ' , 'Dissolved Oxygen Concentration         ',  'mol-C/L' , .true.  , .false., .false., .false.  
    28    sn_tracer(4)   = 'CaCO3   ' , 'Calcite Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. 
    29    sn_tracer(5)   = 'PO4     ' , 'Phosphate Concentration                ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    30    sn_tracer(6)   = 'POC     ' , 'Small organic carbon Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. 
    31    sn_tracer(7)   = 'Si      ' , 'Silicate Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    32    sn_tracer(8)   = 'PHY     ' , 'Nanophytoplankton Concentration        ',  'mol-C/L' , .false. , .false., .false., .false. 
    33    sn_tracer(9)   = 'ZOO     ' , 'Microzooplankton Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. 
    34    sn_tracer(10)  = 'DOC     ' , 'Dissolved organic Concentration        ',  'mol-C/L' , .true.  , .false., .true. , .false. 
    35    sn_tracer(11)  = 'PHY2    ' , 'Diatoms Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. 
    36    sn_tracer(12)  = 'ZOO2    ' , 'Mesozooplankton Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. 
    37    sn_tracer(13)  = 'DSi     ' , 'Diatoms Silicate Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. 
    38    sn_tracer(14)  = 'Fer     ' , 'Dissolved Iron Concentration           ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    39    sn_tracer(15)  = 'BFe     ' , 'Big iron particles Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    40    sn_tracer(16)  = 'GOC     ' , 'Big organic carbon Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    41    sn_tracer(17)  = 'SFe     ' , 'Small iron particles Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. 
    42    sn_tracer(18)  = 'DFe     ' , 'Diatoms iron  Concentration            ',  'mol-C/L' , .false. , .false., .false., .false. 
    43    sn_tracer(19)  = 'GSi     ' , 'Sinking biogenic Silicate Concentration',  'mol-C/L' , .false. , .false., .false., .false. 
    44    sn_tracer(20)  = 'NFe     ' , 'Nano iron Concentration                ',  'mol-C/L' , .false. , .false., .false., .false. 
    45    sn_tracer(21)  = 'NCHL    ' , 'Nano chlorophyl Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. 
    46    sn_tracer(22)  = 'DCHL    ' , 'Diatoms chlorophyl Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    47    sn_tracer(23)  = 'NO3     ' , 'Nitrates Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    48    sn_tracer(24)  = 'NH4     ' , 'Ammonium Concentration                 ',  'mol-C/L' , .false. , .false., .false., .false. 
     24!                !    name   !           title of the field              !   units     ! init    ! sbc    ! cbc    !  obc    !  ais 
     25   sn_tracer(1)   = 'DIC     ' , 'Dissolved inorganic Concentration      ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false. 
     26   sn_tracer(2)   = 'Alkalini' , 'Total Alkalinity Concentration         ',  'eq/L '   , .true.  , .false., .true. , .false. , .false. 
     27   sn_tracer(3)   = 'O2      ' , 'Dissolved Oxygen Concentration         ',  'mol-C/L' , .true.  , .false., .false., .false. , .false. 
     28   sn_tracer(4)   = 'CaCO3   ' , 'Calcite Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     29   sn_tracer(5)   = 'PO4     ' , 'Phosphate Concentration                ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     30   sn_tracer(6)   = 'POC     ' , 'Small organic carbon Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     31   sn_tracer(7)   = 'Si      ' , 'Silicate Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     32   sn_tracer(8)   = 'PHY     ' , 'Nanophytoplankton Concentration        ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     33   sn_tracer(9)   = 'ZOO     ' , 'Microzooplankton Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     34   sn_tracer(10)  = 'DOC     ' , 'Dissolved organic Concentration        ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false. 
     35   sn_tracer(11)  = 'PHY2    ' , 'Diatoms Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     36   sn_tracer(12)  = 'ZOO2    ' , 'Mesozooplankton Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     37   sn_tracer(13)  = 'DSi     ' , 'Diatoms Silicate Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     38   sn_tracer(14)  = 'Fer     ' , 'Dissolved Iron Concentration           ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .true. 
     39   sn_tracer(15)  = 'BFe     ' , 'Big iron particles Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     40   sn_tracer(16)  = 'GOC     ' , 'Big organic carbon Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     41   sn_tracer(17)  = 'SFe     ' , 'Small iron particles Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     42   sn_tracer(18)  = 'DFe     ' , 'Diatoms iron  Concentration            ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     43   sn_tracer(19)  = 'GSi     ' , 'Sinking biogenic Silicate Concentration',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     44   sn_tracer(20)  = 'NFe     ' , 'Nano iron Concentration                ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     45   sn_tracer(21)  = 'NCHL    ' , 'Nano chlorophyl Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     46   sn_tracer(22)  = 'DCHL    ' , 'Diatoms chlorophyl Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     47   sn_tracer(23)  = 'NO3     ' , 'Nitrates Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     48   sn_tracer(24)  = 'NH4     ' , 'Ammonium Concentration                 ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
    4949/ 
    5050!----------------------------------------------------------------------- 
    … …  
    5757!          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    5858!          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
    59    sn_trcdta(1)  = 'data_DIC_nomask'        ,        -12.       ,  'DIC'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    60    sn_trcdta(2)  = 'data_Alkalini_nomask'   ,        -12.       ,  'Alkalini',    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    61    sn_trcdta(3)  = 'data_O2_nomask'         ,        -1.        ,  'O2'      ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    62    sn_trcdta(5)  = 'data_PO4_nomask'        ,        -1.        ,  'PO4'     ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    63    sn_trcdta(7)  = 'data_Si_nomask'         ,        -1.        ,  'Si'      ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    64    sn_trcdta(10) = 'data_DOC_nomask'        ,        -12.       ,  'DOC'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    65    sn_trcdta(14) = 'data_Fer_nomask'        ,        -12.       ,  'Fer'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    66    sn_trcdta(23) = 'data_NO3_nomask'        ,        -1.        ,  'NO3'     ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    67    rn_trfac(1)   =   1.0e-06  !  multiplicative factor 
    68    rn_trfac(2)   =   1.0e-06  !  -      -      -     - 
     59   sn_trcdta(1)  = 'data_DIC_nomask.nc',        -12        ,  'PiDIC'  ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     60   sn_trcdta(2)  = 'data_ALK_nomask.nc',        -12        ,  'TALK'   ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     61   sn_trcdta(3)  = 'data_OXY_nomask.nc',        -1         ,  'O2'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     62   sn_trcdta(5)  = 'data_PO4_nomask.nc',        -1         ,  'PO4'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     63   sn_trcdta(7)  = 'data_SIL_nomask.nc',        -1         ,  'Si'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     64   sn_trcdta(10) = 'data_DOC_nomask.nc',        -1         ,  'DOC'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     65   sn_trcdta(14) = 'data_FER_nomask.nc',        -1         ,  'Fer'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     66   sn_trcdta(23) = 'data_NO3_nomask.nc',        -1         ,  'NO3'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     67   rn_trfac(1)   =   1.028e-06  !  multiplicative factor 
     68   rn_trfac(2)   =   1.028e-06  !  -      -      -     - 
    6969   rn_trfac(3)   =  44.6e-06  !  -      -      -     - 
    7070   rn_trfac(5)   = 122.0e-06  !  -      -      -     - 
    7171   rn_trfac(7)   =   1.0e-06  !  -      -      -     - 
    72    rn_trfac(10)  =   1.0      !  -      -      -     - 
    73    rn_trfac(14)  =   1.0      !  -      -      -     - 
     72   rn_trfac(10)  =   1.0e-06  !  -      -      -     - 
     73   rn_trfac(14)  =   1.0e-06  !  -      -      -     - 
    7474   rn_trfac(23)  =   7.6e-06  !  -      -      -     - 
    7575/ 
    … …  
    110110!                !  file name        ! frequency (hours) ! variable      ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    111111!                !                   !  (if <0  months)  !   name        !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
    112    sn_trcsbc(5)  = 'dust.orca.new'   ,       -1          , 'dustpo4'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    113    sn_trcsbc(7)  = 'dust.orca.new'   ,       -1          , 'dustsi'      ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    114    sn_trcsbc(14) = 'dust.orca.new'   ,       -1          , 'dustfer'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    115    sn_trcsbc(23) = 'ndeposition.orca',      -12          , 'ndep'        ,  .false.     , .true. , 'yearly'  , ''       , ''    , '' 
     112   sn_trcsbc(5)  = 'dust.orca.new'   ,       -1          , 'dustpo4'     ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     113   sn_trcsbc(7)  = 'dust.orca.new'   ,       -1          , 'dustsi'      ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     114   sn_trcsbc(14) = 'dust.orca.new'   ,       -1          , 'dustfer'     ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     115   sn_trcsbc(23) = 'ndeposition.orca',      -12          , 'ndep2'       ,  .false.     , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
    116116   rn_trsfac(5)  = 8.264e-02   !  (  0.021 / 31. * 122 ) 
    117117   rn_trsfac(7)  = 3.313e-01     !  ( 8.8   / 28.1 ) 
    … …  
    120120   rn_sbc_time   =  1.          !  Time scaling factor for SBC and CBC data (seconds in a day) 
    121121   ! 
    122    sn_trccbc(1)  = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    123    sn_trccbc(2)  = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    124    sn_trccbc(5)  = 'river.orca'      ,    120            , 'riverdip'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    125    sn_trccbc(7)  = 'river.orca'      ,    120            , 'riverdsi'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    126    sn_trccbc(10) = 'river.orca'      ,    120            , 'riverdoc'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    127    sn_trccbc(14) = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    128    sn_trccbc(23) = 'river.orca'      ,    120            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     122   sn_trccbc(1)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     123   sn_trccbc(2)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     124   sn_trccbc(5)  = 'river.orca'      ,    -12            , 'riverdip'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     125   sn_trccbc(7)  = 'river.orca'      ,    -12            , 'riverdsi'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     126   sn_trccbc(10) = 'river.orca'      ,    -12            , 'riverdoc'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     127   sn_trccbc(14) = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     128   sn_trccbc(23) = 'river.orca'      ,    -12            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    129129   rn_trcfac(1)  = 8.333e+01   !  ( data in Mg/m2/yr : 1e3/12/ryyss) 
    130130   rn_trcfac(2)  = 8.333e+01   !  ( 1e3 /12 ) 
    … …  
    140140!----------------------------------------------------------------------- 
    141141/ 
     142!----------------------------------------------------------------------- 
     143&namtrc_ais      !  Representation of Antarctic Ice Sheet tracers supply 
     144!----------------------------------------------------------------------- 
     145/ 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg

    r12845 r14054  
    1919   ln_c14        =  .false. 
    2020! 
    21    ln_trcdta     =  .true.   !  Initialisation from data input file (T) or not (F) 
    22    ln_trcbc      =  .false.  !  Enables Boundary conditions 
     21   ln_trcdta     =  .true.  !  Initialisation from data input file (T) or not (F) 
     22   ln_trcbc      =  .true.  !  Enables Boundary conditions 
    2323!                !           !                                           !             !         ! 
    24 !                !    name   !           title of the field              !   units     ! init    ! sbc    ! cbc    !  obc  !  
    25    sn_tracer(1)   = 'DIC     ' , 'Dissolved inorganic Concentration      ',  'mol-C/L' , .true.  , .false., .true. , .false.  
    26    sn_tracer(2)   = 'Alkalini' , 'Total Alkalinity Concentration         ',  'eq/L '   , .true.  , .false., .true. , .false. 
    27    sn_tracer(3)   = 'O2      ' , 'Dissolved Oxygen Concentration         ',  'mol-C/L' , .true.  , .false., .false., .false.  
    28    sn_tracer(4)   = 'CaCO3   ' , 'Calcite Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. 
    29    sn_tracer(5)   = 'PO4     ' , 'Phosphate Concentration                ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    30    sn_tracer(6)   = 'POC     ' , 'Small organic carbon Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. 
    31    sn_tracer(7)   = 'Si      ' , 'Silicate Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    32    sn_tracer(8)   = 'PHY     ' , 'Nanophytoplankton Concentration        ',  'mol-C/L' , .false. , .false., .false., .false. 
    33    sn_tracer(9)   = 'ZOO     ' , 'Microzooplankton Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. 
    34    sn_tracer(10)  = 'DOC     ' , 'Dissolved organic Concentration        ',  'mol-C/L' , .true.  , .false., .true. , .false. 
    35    sn_tracer(11)  = 'PHY2    ' , 'Diatoms Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. 
    36    sn_tracer(12)  = 'ZOO2    ' , 'Mesozooplankton Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. 
    37    sn_tracer(13)  = 'DSi     ' , 'Diatoms Silicate Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. 
    38    sn_tracer(14)  = 'Fer     ' , 'Dissolved Iron Concentration           ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    39    sn_tracer(15)  = 'BFe     ' , 'Big iron particles Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    40    sn_tracer(16)  = 'GOC     ' , 'Big organic carbon Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    41    sn_tracer(17)  = 'SFe     ' , 'Small iron particles Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. 
    42    sn_tracer(18)  = 'DFe     ' , 'Diatoms iron  Concentration            ',  'mol-C/L' , .false. , .false., .false., .false. 
    43    sn_tracer(19)  = 'GSi     ' , 'Sinking biogenic Silicate Concentration',  'mol-C/L' , .false. , .false., .false., .false. 
    44    sn_tracer(20)  = 'NFe     ' , 'Nano iron Concentration                ',  'mol-C/L' , .false. , .false., .false., .false. 
    45    sn_tracer(21)  = 'NCHL    ' , 'Nano chlorophyl Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. 
    46    sn_tracer(22)  = 'DCHL    ' , 'Diatoms chlorophyl Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. 
    47    sn_tracer(23)  = 'NO3     ' , 'Nitrates Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. 
    48    sn_tracer(24)  = 'NH4     ' , 'Ammonium Concentration                 ',  'mol-C/L' , .false. , .false., .false., .false. 
     24!                !    name   !           title of the field              !   units     ! init    ! sbc    ! cbc    !  obc    !  ais 
     25   sn_tracer(1)   = 'DIC     ' , 'Dissolved inorganic Concentration      ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false. 
     26   sn_tracer(2)   = 'Alkalini' , 'Total Alkalinity Concentration         ',  'eq/L '   , .true.  , .false., .true. , .false. , .false. 
     27   sn_tracer(3)   = 'O2      ' , 'Dissolved Oxygen Concentration         ',  'mol-C/L' , .true.  , .false., .false., .false. , .false. 
     28   sn_tracer(4)   = 'CaCO3   ' , 'Calcite Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     29   sn_tracer(5)   = 'PO4     ' , 'Phosphate Concentration                ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     30   sn_tracer(6)   = 'POC     ' , 'Small organic carbon Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     31   sn_tracer(7)   = 'Si      ' , 'Silicate Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     32   sn_tracer(8)   = 'PHY     ' , 'Nanophytoplankton Concentration        ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     33   sn_tracer(9)   = 'ZOO     ' , 'Microzooplankton Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     34   sn_tracer(10)  = 'DOC     ' , 'Dissolved organic Concentration        ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false. 
     35   sn_tracer(11)  = 'PHY2    ' , 'Diatoms Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     36   sn_tracer(12)  = 'ZOO2    ' , 'Mesozooplankton Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     37   sn_tracer(13)  = 'DSi     ' , 'Diatoms Silicate Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     38   sn_tracer(14)  = 'Fer     ' , 'Dissolved Iron Concentration           ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .true. 
     39   sn_tracer(15)  = 'BFe     ' , 'Big iron particles Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     40   sn_tracer(16)  = 'GOC     ' , 'Big organic carbon Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     41   sn_tracer(17)  = 'SFe     ' , 'Small iron particles Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     42   sn_tracer(18)  = 'DFe     ' , 'Diatoms iron  Concentration            ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     43   sn_tracer(19)  = 'GSi     ' , 'Sinking biogenic Silicate Concentration',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     44   sn_tracer(20)  = 'NFe     ' , 'Nano iron Concentration                ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     45   sn_tracer(21)  = 'NCHL    ' , 'Nano chlorophyl Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     46   sn_tracer(22)  = 'DCHL    ' , 'Diatoms chlorophyl Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
     47   sn_tracer(23)  = 'NO3     ' , 'Nitrates Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false. 
     48   sn_tracer(24)  = 'NH4     ' , 'Ammonium Concentration                 ',  'mol-C/L' , .false. , .false., .false., .false. , .false. 
    4949/ 
    5050!----------------------------------------------------------------------- 
    … …  
    5757!          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    5858!          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
    59    sn_trcdta(1)  = 'data_DIC_nomask'        ,        -12.       ,  'DIC'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    60    sn_trcdta(2)  = 'data_Alkalini_nomask'   ,        -12.       ,  'Alkalini',    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    61    sn_trcdta(3)  = 'data_O2_nomask'         ,        -1.        ,  'O2'      ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    62    sn_trcdta(5)  = 'data_PO4_nomask'        ,        -1.        ,  'PO4'     ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    63    sn_trcdta(7)  = 'data_Si_nomask'         ,        -1.        ,  'Si'      ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    64    sn_trcdta(10) = 'data_DOC_nomask'        ,        -12.       ,  'DOC'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    65    sn_trcdta(14) = 'data_Fer_nomask'        ,        -12.       ,  'Fer'     ,    .false.   , .true. , 'yearly'  , ''       , ''   , '' 
    66    sn_trcdta(23) = 'data_NO3_nomask'        ,        -1.        ,  'NO3'     ,    .true.    , .true. , 'yearly'  , ''       , ''   , '' 
    67    rn_trfac(1)   =   1.0e-06  !  multiplicative factor 
    68    rn_trfac(2)   =   1.0e-06  !  -      -      -     - 
     59   sn_trcdta(1)  = 'data_DIC_nomask.nc',        -12        ,  'PiDIC'  ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     60   sn_trcdta(2)  = 'data_ALK_nomask.nc',        -12        ,  'TALK'   ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     61   sn_trcdta(3)  = 'data_OXY_nomask.nc',        -1         ,  'O2'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     62   sn_trcdta(5)  = 'data_PO4_nomask.nc',        -1         ,  'PO4'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     63   sn_trcdta(7)  = 'data_SIL_nomask.nc',        -1         ,  'Si'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     64   sn_trcdta(10) = 'data_DOC_nomask.nc',        -1         ,  'DOC'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     65   sn_trcdta(14) = 'data_FER_nomask.nc',        -1         ,  'Fer'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     66   sn_trcdta(23) = 'data_NO3_nomask.nc',        -1         ,  'NO3'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , '' 
     67   rn_trfac(1)   =   1.028e-06  !  multiplicative factor 
     68   rn_trfac(2)   =   1.028e-06  !  -      -      -     - 
    6969   rn_trfac(3)   =  44.6e-06  !  -      -      -     - 
    7070   rn_trfac(5)   = 122.0e-06  !  -      -      -     - 
    7171   rn_trfac(7)   =   1.0e-06  !  -      -      -     - 
    72    rn_trfac(10)  =   1.0      !  -      -      -     - 
    73    rn_trfac(14)  =   1.0      !  -      -      -     - 
     72   rn_trfac(10)  =   1.0e-06  !  -      -      -     - 
     73   rn_trfac(14)  =   1.0e-06  !  -      -      -     - 
    7474   rn_trfac(23)  =   7.6e-06  !  -      -      -     - 
    7575/ 
    … …  
    110110!                !  file name        ! frequency (hours) ! variable      ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    111111!                !                   !  (if <0  months)  !   name        !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
    112    sn_trcsbc(5)  = 'dust.orca.new'   ,       -1          , 'dustpo4'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    113    sn_trcsbc(7)  = 'dust.orca.new'   ,       -1          , 'dustsi'      ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    114    sn_trcsbc(14) = 'dust.orca.new'   ,       -1          , 'dustfer'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    115    sn_trcsbc(23) = 'ndeposition.orca',      -12          , 'ndep'        ,  .false.     , .true. , 'yearly'  , ''       , ''    , '' 
     112   sn_trcsbc(5)  = 'dust.orca.new'   ,       -1          , 'dustpo4'     ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     113   sn_trcsbc(7)  = 'dust.orca.new'   ,       -1          , 'dustsi'      ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     114   sn_trcsbc(14) = 'dust.orca.new'   ,       -1          , 'dustfer'     ,  .true.      , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
     115   sn_trcsbc(23) = 'ndeposition.orca',      -12          , 'ndep2'       ,  .false.     , .true. , 'yearly'  , 'weights_2D_r360x180_bilin.nc'       , ''    , '' 
    116116   rn_trsfac(5)  = 8.264e-02   !  (  0.021 / 31. * 122 ) 
    117117   rn_trsfac(7)  = 3.313e-01     !  ( 8.8   / 28.1 ) 
    … …  
    120120   rn_sbc_time   =  1.          !  Time scaling factor for SBC and CBC data (seconds in a day) 
    121121   ! 
    122    sn_trccbc(1)  = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    123    sn_trccbc(2)  = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    124    sn_trccbc(5)  = 'river.orca'      ,    120            , 'riverdip'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    125    sn_trccbc(7)  = 'river.orca'      ,    120            , 'riverdsi'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    126    sn_trccbc(10) = 'river.orca'      ,    120            , 'riverdoc'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    127    sn_trccbc(14) = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    128    sn_trccbc(23) = 'river.orca'      ,    120            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     122   sn_trccbc(1)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     123   sn_trccbc(2)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     124   sn_trccbc(5)  = 'river.orca'      ,    -12            , 'riverdip'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     125   sn_trccbc(7)  = 'river.orca'      ,    -12            , 'riverdsi'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     126   sn_trccbc(10) = 'river.orca'      ,    -12            , 'riverdoc'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     127   sn_trccbc(14) = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
     128   sn_trccbc(23) = 'river.orca'      ,    -12            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , '' 
    129129   rn_trcfac(1)  = 8.333e+01   !  ( data in Mg/m2/yr : 1e3/12/ryyss) 
    130130   rn_trcfac(2)  = 8.333e+01   !  ( 1e3 /12 ) 
    … …  
    140140!----------------------------------------------------------------------- 
    141141/ 
     142!----------------------------------------------------------------------- 
     143&namtrc_ais      !  Representation of Antarctic Ice Sheet tracers supply 
     144!----------------------------------------------------------------------- 
     145/ 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/axis_def_nemo.xml

    r12377 r14054  
    1414      <axis id="depthv"  long_name="Vertical V levels" unit="m" positive="down" /> 
    1515      <axis id="depthw"  long_name="Vertical W levels" unit="m" positive="down" /> 
     16      <axis id="depthf"  long_name="Vertical F levels" unit="m" positive="down" /> 
    1617      <axis id="nfloat"  long_name="Float number"      unit="-"                 /> 
    1718      <axis id="icbcla"  long_name="Iceberg class"      unit="1"               /> 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/domain_def_nemo.xml

    r12276 r14054  
    181181     <domain id="EqW" domain_ref="grid_W" > <zoom_domain id="EqW"/> </domain> 
    182182 
     183     <!--   F grid   --> 
     184     <domain id="grid_F" long_name="grid F"/> 
     185      
    183186              <!--   zonal mean grid   --> 
    184187     <domain_group id="gznl"> 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/field_def_nemo-oce.xml

    r14049 r14054  
    171171        <field id="tosmint_pot"  long_name="vertical integral of potential temperature times density"   standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" /> 
    172172 
    173  
     173        <field id="ht"           long_name="water column height at T point"                     standard_name="water_column_height_T"                      unit="m" /> 
    174174        <field id="ssh"          long_name="sea surface height"                                 standard_name="sea_surface_height_above_geoid"             unit="m" /> 
    175175        <field id="ssh2"         long_name="square of sea surface height"                       standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field > 
    … …  
    190190 
    191191        <!-- Energy - horizontal divergence --> 
     192        <field id="sKE"          long_name="surface kinetic energy"  standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_2D" /> 
    192193        <field id="hdiv"         long_name="horizontal divergence"                                                          unit="s-1"    grid_ref="grid_T_3D" /> 
    193194 
    … …  
    270271 
    271272      <field_group id="OSMOSIS_T" grid_ref="grid_T_2D"> 
     273        <field id="hml"                 long_name="mixed layr depth"                         unit="m"       /> 
     274        <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
     275        <field id="dh"                  long_name="Pycnocline thickness"                     unit=" m"      /> 
     276        <field id="ibld"                long_name="index of boundary layer depth"            unit="#"       /> 
     277        <field id="imld"                long_name="index of mixed layer depth"            unit="#"       /> 
     278        <field id="zhbl"                long_name="boundary layer depth -grid"                     unit="m"       /> 
     279        <field id="zhml"                long_name="mixed layer depth - grid"                        unit="m"       /> 
     280        <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
     281        <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
     282        <field id="us_x"        long_name="i component of active Stokes drift"                      unit="m/s"     /> 
     283        <field id="us_y"        long_name="j component of active Stokes drift"                      unit="m/s"     /> 
     284        <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    272285        <field id="zwth0"               long_name="surface non-local temperature flux"       unit="deg m/s" /> 
    273286        <field id="zws0"                long_name="surface non-local salinity flux"          unit="psu m/s" /> 
    274         <field id="hbl"                 long_name="boundary layer depth"                     unit="m"       /> 
    275         <field id="hbli"                long_name="initial boundary layer depth"             unit="m"       /> 
    276         <field id="dstokes"             long_name="stokes drift  depth scale"                unit="m"       /> 
    277         <field id="zustke"              long_name="magnitude of stokes drift  at T-points"   unit="m/s"     /> 
    278287        <field id="zwstrc"              long_name="convective velocity scale"                unit="m/s"     /> 
     288        <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    279289        <field id="zwstrl"              long_name="langmuir velocity scale"                  unit="m/s"     /> 
    280         <field id="zustar"              long_name="friction velocity"                        unit="m/s"     /> 
    281         <field id="zhbl"                long_name="boundary layer depth"                     unit="m"       /> 
    282         <field id="zhml"                long_name="mixed layer depth"                        unit="m"       /> 
     290        <field id="zvstr"               long_name="mixed velocity scale"                     unit="m/s"     /> 
     291        <field id="zla"                 long_name="langmuir number"                          unit="m/s"     /> 
    283292        <field id="wind_wave_abs_power" long_name="\rho |U_s| x u*^2"                        unit="mW"      /> 
    284293        <field id="wind_wave_power"     long_name="U_s \dot  tau"                            unit="mW"      /> 
    285294        <field id="wind_power"          long_name="\rho  u*^3"                               unit="mW"      /> 
    286295 
    287         <!-- extra OSMOSIS diagnostics --> 
     296       <!-- interior BL OSMOSIS diagnostics --> 
    288297        <field id="zwthav"              long_name="av turb flux of T in ml"                  unit="deg m/s" /> 
    289298        <field id="zt_ml"               long_name="av T in ml"                               unit="deg"     /> 
     299        <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
     300        <field id="zws_ent"            long_name="entrainment turb flux of S"                unit="10^-3 m/s" /> 
    290301        <field id="zwth_ent"            long_name="entrainment turb flux of T"               unit="deg m/s" /> 
    291         <field id="zhol"                long_name="Hoenekker number"                         unit="#"       /> 
    292         <field id="zdh"                 long_name="Pycnocline  depth - grid"                 unit=" m"      /> 
    293       </field_group> 
    294  
    295       <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" operation="instant" > 
     302        <field id="zwb_ent"            long_name="entrainment turb flux of buoyancy"         unit="m^2/s^-3" /> 
     303  
     304        <field id="zdt_bl"             long_name="temperature jump at base of BL"                 unit="deg"      /> 
     305        <field id="zds_bl"             long_name="salinity jump at base of BL"                 unit="10^-3"      /> 
     306        <field id="zdb_bl"             long_name="buoyancy jump at base of BL"                 unit="m/s^2"      /> 
     307        <field id="zdu_bl"             long_name="u jump at base of BL"                       unit="m/s"      /> 
     308        <field id="zdv_bl"             long_name="v jump at base of BL"                       unit="m/s"      /> 
     309 
     310        <!-- extra OSMOSIS diagnostics for debugging --> 
     311       <field id="zsc_uw_1_0"       long_name="zsc u-momentum flux on T after Stokes"                       unit="m^2/s^2" /> 
     312        <field id="zsc_uw_1_f"       long_name="zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     313        <field id="zsc_vw_1_f"       long_name="zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     314        <field id="zsc_uw_2_f"       long_name="2nd zsc u-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     315        <field id="zsc_vw_2_f"       long_name="2nd zsc v-momentum flux on T after Coriolis"                       unit="m^2/s^2" /> 
     316        <field id="zuw_bse"       long_name="base u-flux T-points"                          unit="m^2/s^2" /> 
     317        <field id="zvw_bse"       long_name="base v-flux T-points"                          unit="m^2/s^2" /> 
     318 
     319       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     320         <field id="hmle"          long_name="OBL FK-layer thickness"                                     unit="m"        /> 
     321        <field id="mld_prof"              long_name="FK-layer depth index"                  unit="#" /> 
     322        <field id="zmld"          long_name="target FK-layer thickness"                                     unit="m"        /> 
     323        <field id="zwb_fk"          long_name="FK b-flux"                                     unit="m^2 s^-3"        /> 
     324        <field id="zwb_fk_b"          long_name="layer averaged FK b-flux"                 unit="m^2 s^-3"       /> 
     325        <field id="zdiff_mle"          long_name="max FK diffusivity in MLE"       unit=" 10^-4 m^2 s^-1"       /> 
     326        <field id="zvel_mle"          long_name="FK velocity scale in MLE"       unit=" m s^-1"       /> 
     327    </field_group> 
     328 
     329      <field_group id="OSMOSIS_W" grid_ref="grid_W_3D" > 
     330        <field id="zviscos"       long_name="BL viscosity"   unit="m^2/s" /> 
    296331        <field id="ghamt"       long_name="non-local temperature flux"                       unit="deg m/s" /> 
    297332        <field id="ghams"       long_name="non-local salinity flux"                          unit="psu m/s" /> 
    298333        <field id="zdtdz_pyc"   long_name="Pycnocline temperature gradient"                  unit=" deg/m"  /> 
    299       </field_group> 
     334        <field id="zdsdz_pyc"   long_name="Pycnocline salinity gradient"                  unit=" 10^-3/m"  /> 
     335        <field id="zdbdz_pyc"   long_name="Pycnocline buoyancy gradient"                  unit=" s^-2"  /> 
     336        <field id="zdudz_pyc"   long_name="Pycnocline u gradient"                  unit=" s^-2"  /> 
     337        <field id="zdvdz_pyc"   long_name="Pycnocline v gradient"                  unit=" s^-2"  /> 
     338 
     339        <!-- extra OSMOSIS diagnostics for debugging --> 
     340         <field id="ghamu_00"       long_name="initial non-local u-momentum flux"   unit="m^2/s^2" /> 
     341        <field id="ghamv_00"       long_name="initial non-local v-momentum flux"   unit="m^2/s^2" /> 
     342        <field id="ghamu_0"       long_name="after dstokes non-local u-momentum flux"   unit="m^2/s^2" /> 
     343        <field id="ghamu_f"       long_name="after Coriolis non-local u-momentum flux"   unit="m^2/s^2" /> 
     344        <field id="ghamv_f"       long_name="after Coriolis  non-local v-momentum flux"   unit="m^2/s^2" /> 
     345        <field id="ghamu_b"       long_name="after buoyancy added non-local u-momentum flux"   unit="m^2/s^2" /> 
     346        <field id="ghamv_b"       long_name="after buoyancy added  non-local v-momentum flux"  unit="m^2/s^2" /> 
     347        <field id="ghamu_1"       long_name="after entrainment non-local u-momentum flux"   unit="m^2/s^2" /> 
     348        <field id="ghamv_1"       long_name="after entrainment  non-local v-momentum flux"  unit="m^2/s^2" /> 
     349     </field_group> 
    300350 
    301351      <field_group id="OSMOSIS_U" grid_ref="grid_U_2D" > 
    302352        <field id="ghamu"       long_name="non-local u-momentum flux"   grid_ref="grid_U_3D" unit="m^2/s^2" /> 
    303         <field id="us_x"        long_name="i component of Stokes drift"                      unit="m/s"     /> 
    304       </field_group> 
     353       <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     354       <field id="zdtdx"          long_name="FK  T x-gradient"                                     unit=" deg C m^-1"        /> 
     355        <field id="zdsdx"          long_name="FK  S x-gradient"                                     unit=" 10^-3 m^-1"        /> 
     356        <field id="dbdx_mle"          long_name="FK  B x-gradient"                                     unit=" s^-2"        /> 
     357     </field_group> 
    305358 
    306359      <field_group id="OSMOSIS_V" grid_ref="grid_V_2D" > 
    307360        <field id="ghamv"       long_name="non-local v-momentum flux"   grid_ref="grid_V_3D" unit="m^2/s^2" /> 
    308         <field id="us_y"        long_name="j component of Stokes drift"                      unit="m/s"     /> 
     361        <!-- FK_OSM OSMOSIS diagnostics (require also ln_osm_mle=.true.--> 
     362        <field id="zdtdy"          long_name="FK T y-gradient"                                     unit=" deg C m^-1"        /> 
     363        <field id="zdsdy"          long_name="FK S y-gradient"                                     unit=" 10^-3 m^-1"        /> 
     364        <field id="dbdy_mle"          long_name="FK B y-gradient"                                     unit=" s^-2"        /> 
    309365      </field_group> 
    310366 
    … …  
    501557 
    502558      <field_group id="grid_U"   grid_ref="grid_U_2D"> 
     559        <field id="hu"            long_name="water column height at U point"                         standard_name="water_column_height_U"       unit="m" /> 
    503560        <field id="e2u"           long_name="U-cell width in meridional direction"                   standard_name="cell_width"                  unit="m"                               /> 
    504561        <field id="e3u"           long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" /> 
    … …  
    571628        <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" /> 
    572629        <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D" /> 
     630        <field id="hv"            long_name="water column height at V point"                        standard_name="water_column_height_V"       unit="m" /> 
    573631        <field id="vtau"         long_name="Wind Stress along j-axis"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            /> 
    574632        <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
    … …  
    679737 
    680738      <!-- F grid --> 
     739      <field_group id="grid_F" grid_ref="grid_F_2D"> 
     740   <field id="e3f"          long_name="F-cell thickness"                    standard_name="cell_thickness"        unit="m"   grid_ref="grid_F_3D" /> 
     741   <field id="e3f_0"        long_name="F-cell thickness"                    standard_name="cell_thickness"        unit="m"   grid_ref="grid_F_3D" /> 
     742        <field id="hf"           long_name="water column height at F point"    standard_name="water_column_height_F"  unit="m"                     /> 
     743        <field id="sKEf"         long_name="surface kinetic energy at F point" standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2" /> 
     744        <field id="relvor"       long_name="relative vorticity"                standard_name="relative_vorticity"     unit="1/s"                   /> 
     745        <field id="plavor"       long_name="planetary vorticity"               standard_name="planetary_vorticity"    unit="1/s"                   /> 
     746        <field id="relpotvor"    long_name="relative potential vorticity"      standard_name="relpot_vorticity"       unit="1/m.s"                 /> 
     747        <field id="abspotvor"    long_name="absolute potential vorticity"      standard_name="abspot_vorticity"       unit="1/m.s"                 /> 
     748        <field id="Ens"          long_name="enstrophy"                         standard_name="enstrophy"              unit="1/m2.s2"               /> 
     749      </field_group>  
     750  
    681751      <!-- AGRIF sponge --> 
    682752      <field id="agrif_spf"    long_name=" AGRIF f-sponge coefficient"   unit=" " /> 
    … …  
    841911     <field id="strd_zdfp"     long_name="salinity   -trend: pure vert. diffusion"   unit="1e-3/s" /> 
    842912 
    843      <!-- --> 
     913     <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     914     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s" /> 
     915     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s" /> 
     916 
     917 
     918    <!-- --> 
    844919     <field id="ttrd_dmp"      long_name="temperature-trend: interior restoring"        unit="degC/s" /> 
    845920     <field id="strd_dmp"      long_name="salinity   -trend: interior restoring"        unit="1e-3/s" /> 
    … …  
    877952     <field id="strd_zdfp_e3t"     unit="1e-3/s * m"  >  strd_zdfp * e3t </field> 
    878953 
     954          <!-- ln_zdfosm=T only (OSMOSIS-OBL) --> 
     955     <field id="ttrd_osm_e3t"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field> 
     956     <field id="strd_osm_e3t"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field> 
     957      
    879958     <!-- --> 
    880959     <field id="ttrd_dmp_e3t"      unit="degC/s * m"  >  ttrd_dmp * e3t </field> 
    … …  
    892971     <field id="ttrd_totad_li"    long_name="layer integrated heat-trend: total advection"         unit="W/m^2"     > ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field> 
    893972     <field id="strd_totad_li"    long_name="layer integrated salt-trend: total advection"         unit="kg/(m^2 s)"    > strd_totad_e3t * 1026.0 * 0.001  </field> 
     973     <field id="ttrd_osm_li"    long_name="layer integrated heat-trend: non-local OSM"         unit="W/m^2"     > ttrd_osm_e3t * 1026.0 * 3991.86795711963 </field> 
     974     <field id="strd_osm_li"    long_name="layer integrated salt-trend: non-local OSM"         unit="kg/(m^2 s)"    > strd_osm_e3t * 1026.0 * 0.001  </field> 
    894975     <field id="ttrd_evd_li"      long_name="layer integrated heat-trend: EVD convection"          unit="W/m^2"    > ttrd_evd_e3t * 1026.0 * 3991.86795711963 </field> 
    895976     <field id="strd_evd_li"      long_name="layer integrated salt-trend: EVD convection"          unit="kg/(m^2 s)"  > strd_evd_e3t * 1026.0 * 0.001  </field> 
    … …  
    10991180    </field_group> 
    11001181 
     1182    <!-- TMB diagnostic output --> 
     1183    <field_group  id="1h_grid_T_tmb" grid_ref="grid_T_2D" operation="instant"> 
     1184      <field id="top_temp"           name="votemper_top"  unit="degC"  /> 
     1185      <field id="mid_temp"           name="votemper_mid"  unit="degC"  /> 
     1186      <field id="bot_temp"           name="votemper_bot"  unit="degC"  /> 
     1187      <field id="top_sal"            name="vosaline_top"  unit="psu"   /> 
     1188      <field id="mid_sal"            name="vosaline_mid"  unit="psu"   /> 
     1189      <field id="bot_sal"            name="vosaline_bot"  unit="psu"   /> 
     1190      <field id="sshnmasked"         name="sossheig"      unit="m"     />  
     1191    </field_group> 
     1192 
    11011193    <field_group  id="1h_grid_U_tmb" grid_ref="grid_U_2D" operation="instant"> 
    11021194      <field id="top_u"           name="vozocrtx_top"  unit="m/s"  /> 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/grid_def_nemo.xml

    r12377 r14054  
    5353         <domain domain_ref="grid_W" /> 
    5454         <axis axis_ref="depthw" /> 
     55       </grid> 
     56       <!--  --> 
     57       <grid id="grid_F_2D" > 
     58         <domain domain_ref="grid_F" /> 
     59       </grid> 
     60        <!--  --> 
     61       <grid id="grid_F_3D" > 
     62         <domain domain_ref="grid_F" /> 
     63         <axis axis_ref="depthf" /> 
    5564       </grid> 
    5665        <!--  --> 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/namelist_ref

    r14049 r14054  
    629629   ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0 
    630630   rn_speed_limit          = 0.      ! CFL speed limit for a berg 
    631  
     631   ! 
     632   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 
     633      ln_icb_grd           = .false. ! ground icb when icb bottom level hit oce bottom level (need ln_M2016 to be activated) 
     634   ! 
    632635   cn_dir      = './'      !  root directory for the calving data location 
    633636   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
    … …  
    995998   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) 
    996999   ln_dynvor_een = .false. !  energy & enstrophy scheme 
    997       nn_een_e3f = 0          ! =0  e3f = mi(mj(e3t))/4 
    998       !                       ! =1  e3f = mi(mj(e3t))/mi(mj( tmask)) 
     1000   ! 
    9991001   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T)        ==>>> PLEASE DO NOT ACTIVATE 
    1000       !                    !  (f-point vorticity schemes only) 
     1002   !                       !  (f-point vorticity schemes only) 
     1003   ! 
     1004   nn_e3f_typ = 0          !  type of e3f (EEN, ENE, ENS, MIX only)  =0  e3f = mi(mj(e3t))/4 
     1005   !                       !                                         =1  e3f = mi(mj(e3t))/mi(mj( tmask)) 
    10011006/ 
    10021007!----------------------------------------------------------------------- 
    … …  
    10301035   !                       !  Type of the operator : 
    10311036   ln_dynldf_OFF = .false.     !  No operator (i.e. no explicit diffusion) 
     1037   nn_dynldf_typ = 0           !  =0 div-rot (default)   ;   =1 symmetric 
    10321038   ln_dynldf_lap = .false.     !    laplacian operator 
    10331039   ln_dynldf_blp = .false.     !  bilaplacian operator 
    … …  
    11601166   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F) 
    11611167      nn_mxlice    = 2        ! type of scaling under sea-ice 
    1162                               !    = 0 no scaling under sea-ice 
    1163                               !    = 1 scaling with constant sea-ice thickness 
    1164                               !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
    1165                               !    = 3 scaling with maximum sea-ice thickness 
     1168      !                       !    = 0 no scaling under sea-ice 
     1169      !                       !    = 1 scaling with constant sea-ice thickness 
     1170      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
     1171      !                       !    = 3 scaling with maximum sea-ice thickness 
    11661172      rn_mxlice   = 10.       ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 
    11671173   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value 
    … …  
    11701176      rn_lc       =   0.15    !  coef. associated to Langmuir cells 
    11711177   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to NIWs 
    1172                               !        = 0 none ; = 1 add a tke source below the ML 
    1173                               !        = 2 add a tke source just at the base of the ML 
    1174                               !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T) 
     1178   !                          !        = 0 none ; = 1 add a tke source below the ML 
     1179   !                          !        = 2 add a tke source just at the base of the ML 
     1180   !                          !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T) 
    11751181      rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2) 
    11761182      nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML 
    1177                               !        = 0  constant 10 m length scale 
    1178                               !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
     1183      !                       !        = 0  constant 10 m length scale 
     1184      !                       !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
    11791185   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice 
    11801186   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking 
    … …  
    12131219&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T) 
    12141220!----------------------------------------------------------------------- 
    1215    ln_use_osm_la = .false.      !  Use namelist  rn_osm_la 
     1221   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    12161222   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    1217    rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     1223   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor 
     1224   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4 
     1225   rn_osm_bl_thresh   = 5.e-5      !Threshold buoyancy for deepening of OSBL base 
    12181226   nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv 
    12191227   ln_dia_osm = .true.         ! output OSMOSIS-OBL variables 
    … …  
    12231231   rn_difri  =  0.005          ! max Ri# diffusivity at Ri_g = 0 (m^2/s) 
    12241232   ln_convmix  = .true.        ! Use convective instability mixing below BL 
    1225    rn_difconv = 1.             ! diffusivity when unstable below BL  (m2/s) 
     1233   rn_difconv = 1. !0.01 !1.             ! diffusivity when unstable below BL  (m2/s) 
     1234   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
    12261235   nn_osm_wave = 0             ! Method used to calculate Stokes drift 
    12271236      !                        !  = 2: Use ECMWF wave fields 
    12281237      !                        !  = 1: Pierson Moskowitz wave spectrum 
    12291238      !                        !  = 0: Constant La# = 0.3 
    1230 / 
     1239   nn_osm_SD_reduce = 0        ! Method used to get active Stokes drift from surface value 
     1240      !                        !  = 0: No reduction 
     1241                               !  = 1: use SD avged over top 10% hbl 
     1242                               !  = 2:use surface value of SD fit to slope at rn_osm_hblfrac*hbl below surface 
     1243   ln_zdfosm_ice_shelter = .true.  ! reduce surface SD and depth scale under ice 
     1244   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
     1245/ 
     1246!----------------------------------------------------------------------- 
     1247&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF) 
     1248!----------------------------------------------------------------------- 
     1249   rn_osm_mle_ce       = 0.06      ! magnitude of the MLE (typical value: 0.06 to 0.08) 
     1250   nn_osm_mle          = 0         ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 
     1251   rn_osm_mle_lf       = 5.e+3     ! typical scale of mixed layer front (meters)                      (case rn_osm_mle=0) 
     1252   rn_osm_mle_time     = 172800.   ! time scale for mixing momentum across the mixed layer (seconds)  (case rn_osm_mle=0) 
     1253   rn_osm_mle_lat      = 20.       ! reference latitude (degrees) of MLE coef.                        (case rn_mle=1) 
     1254   rn_osm_mle_rho_c =    0.01      ! delta rho criterion used to calculate MLD for FK 
     1255   rn_osm_mle_thresh  = 0.0005     ! delta b criterion used for FK MLE criterion 
     1256   rn_osm_mle_tau     = 172800.    ! time scale for FK-OSM (seconds)  (case rn_osm_mle=0) 
     1257   ln_osm_hmle_limit   = .false.   ! limit hmle to rn_osm_hmle_limit*hbl 
     1258   rn_osm_hmle_limit   = 1.2 
     1259   / 
    12311260!----------------------------------------------------------------------- 
    12321261&namzdf_mfc     !   Mass Flux Convection 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/namelist_top_ref

    r12377 r14054  
    4141   ln_trcdmp_clo =  .false.  !  damping term (T) or not (F) on closed seas 
    4242   ln_trcbc      =  .false.  !  Surface, Lateral or Open Boundaries conditions 
     43   ln_trcais     =  .false.  !  Antarctic Ice Sheet nutrient supply 
    4344   ! 
    4445   jp_dia3d      = 0         ! Number of 3D diagnostic variables 
    … …  
    149150                             !  = 2 Damping applied to all tracers 
    150151/ 
     152!----------------------------------------------------------------------- 
     153&namtrc_ais      !  Representation of Antarctic Ice Sheet tracers supply 
     154!----------------------------------------------------------------------- 
     155   nn_ais_tr     =  1        !  tracer concentration in iceberg and ice shelf 
     156                             !    = 0 (null concentrations) 
     157                             !    = 1 prescribed concentrations 
     158   rn_icbdep     =  120.     ! Mean underwater depth of iceberg (m) 
     159/ 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/ICE/iceistate.F90

    r14049 r14054  
    2121   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 
    2222   USE eosbn2         ! equation of state 
     23# if defined key_qco 
     24   USE domqco         ! Variable volume 
     25# else 
    2326   USE domvvl         ! Variable volume 
     27# endif 
    2428   USE ice            ! sea-ice: variables 
    2529   USE ice1D          ! sea-ice: thermodynamics variables 
    … …  
    434438         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0 
    435439         ! 
     440#if defined key_qco 
     441         IF( .NOT.ln_linssh )   CALL dom_qco_zgr( Kbb, Kmm )        ! interpolation scale factor, depth and water column 
     442#else 
    436443         IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column 
    437 ! !!st 
    438 !          IF( .NOT.ln_linssh ) THEN 
    439 !             ! 
    440 !             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 
    441 !             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
    442 !             ! 
    443 !             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    444 !                e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 
    445 !                e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    446 !                e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 
    447 !             END DO 
    448 !             ! 
    449 !             ! Reconstruction of all vertical scale factors at now and before time-steps 
    450 !             ! ========================================================================= 
    451 !             ! Horizontal scale factor interpolations 
    452 !             ! -------------------------------------- 
    453 !             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 
    454 !             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 
    455 !             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    456 !             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    457 !             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 
    458 !             ! Vertical scale factor interpolations 
    459 !             ! ------------------------------------ 
    460 !             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
    461 !             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
    462 !             CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
    463 !             CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
    464 !             CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    465 !             ! t- and w- points depth 
    466 !             ! ---------------------- 
    467 !             !!gm not sure of that.... 
    468 !             gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
    469 !             gdepw(:,:,1,Kmm) = 0.0_wp 
    470 !             gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    471 !             DO jk = 2, jpk 
    472 !                gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm) 
    473 !                gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 
    474 !                gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm) 
    475 !             END DO 
    476 !          ENDIF 
     444#endif 
     445 
    477446      ENDIF 
    478447 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/ICE/icerst.F90

    r14049 r14054  
    9999               CALL iom_swap( cxios_context ) 
    100100#else 
    101                clinfo = 'Can not use XIOS in rst_opn' 
    102                CALL ctl_stop(TRIM(clinfo)) 
     101               CALL ctl_stop( 'Can not use XIOS in rst_opn' ) 
    103102#endif 
    104103            ENDIF 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/NST/agrif_oce_interp.F90

    r13286 r14054  
    2828   USE agrif_oce 
    2929   USE phycst 
    30    USE dynspg_ts, ONLY: un_adv, vn_adv 
     30!!!   USE dynspg_ts, ONLY: un_adv, vn_adv 
    3131   ! 
    3232   USE in_out_manager 
    … …  
    5050   INTEGER ::   bdy_tinterp = 0 
    5151 
    52    !!---------------------------------------------------------------------- 
     52   !! * Substitutions 
     53#  include "domzgr_substitute.h90" 
    5354   !! NEMO/NST 4.0 , NEMO Consortium (2018) 
    5455   !! $Id$ 
    … …  
    11921193      !!----------------------------------------------------------------------   
    11931194      IF( before ) THEN 
    1194          IF ( ln_bt_fw ) THEN 
     1195!         IF ( ln_bt_fw ) THEN 
    11951196            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 
    1196          ELSE 
    1197             ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 
    1198          ENDIF 
     1197!         ELSE 
     1198!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 
     1199!         ENDIF 
    11991200      ELSE 
    12001201         zrhot = Agrif_rhot() 
    … …  
    12281229      ! 
    12291230      IF( before ) THEN 
    1230          IF ( ln_bt_fw ) THEN 
     1231!         IF ( ln_bt_fw ) THEN 
    12311232            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
    1232          ELSE 
    1233             ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
    1234          ENDIF 
     1233!         ELSE 
     1234!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1235!         ENDIF 
    12351236      ELSE       
    12361237         zrhot = Agrif_rhot() 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/NST/agrif_oce_sponge.F90

    r13312 r14054  
    3232 
    3333   !! * Substitutions 
     34#  include "domzgr_substitute.h90" 
    3435#  include "do_loop_substitute.h90" 
    3536   !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/NST/agrif_oce_update.F90

    r13286 r14054  
    2727   USE vremap         ! Vertical remapping 
    2828   USE lbclnk  
    29  
     29#if defined key_qco 
     30   USE domqco 
     31#endif 
    3032   IMPLICIT NONE 
    3133   PRIVATE 
    … …  
    3436   PUBLIC   Update_Scales 
    3537 
     38   !! * Substitutions 
     39#  include "domzgr_substitute.h90" 
    3640   !!---------------------------------------------------------------------- 
    3741   !! NEMO/NST 4.0 , NEMO Consortium (2018) 
    … …  
    191195   END SUBROUTINE Agrif_Update_Tke 
    192196 
    193  
    194197   SUBROUTINE Agrif_Update_vvl( ) 
    195198      !!--------------------------------------------- 
    … …  
    201204      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
    202205      ! 
     206#if ! defined key_qco 
    203207      Agrif_UseSpecialValueInUpdate = .TRUE. 
    204208      Agrif_SpecialValueFineGrid = 0. 
    … …  
    213217      CALL dom_vvl_update_UVF 
    214218      CALL Agrif_ParentGrid_To_ChildGrid() 
     219#else 
     220      CALL Agrif_ChildGrid_To_ParentGrid() 
     221      CALL Agrif_Update_qco 
     222      CALL Agrif_ParentGrid_To_ChildGrid() 
     223#endif 
    215224      ! 
    216225   END SUBROUTINE Agrif_Update_vvl 
    217226 
     227 
     228#if defined key_qco 
     229   SUBROUTINE Agrif_Update_qco 
     230      !!--------------------------------------------- 
     231      !!       *** ROUTINE dom_Update_qco *** 
     232      !!--------------------------------------------- 
     233      ! 
     234      ! Save arrays prior update (needed for asselin correction) 
     235      r3t(:,:,Krhs_a) = r3t(:,:,Kmm_a) 
     236      r3u(:,:,Krhs_a) = r3u(:,:,Kmm_a) 
     237      r3v(:,:,Krhs_a) = r3v(:,:,Kmm_a) 
     238 
     239      ! Update r3x arrays from updated ssh 
     240      CALL dom_qco_zgr( Kbb_a, Kmm_a ) 
     241      ! 
     242   END SUBROUTINE Agrif_Update_qco 
     243#endif 
     244 
     245 
     246#if ! defined key_qco 
    218247   SUBROUTINE dom_vvl_update_UVF 
    219248      !!--------------------------------------------- 
    … …  
    224253      REAL(wp):: zcoef 
    225254      !!--------------------------------------------- 
    226  
    227255      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 
    228256                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
    … …  
    290318      ! 
    291319   END SUBROUTINE dom_vvl_update_UVF 
     320#endif 
    292321 
    293322#if defined key_vertical 
    … …  
    13321361   END SUBROUTINE updateAVM 
    13331362 
     1363#if ! defined key_qco 
    13341364   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) 
    13351365      !!--------------------------------------------- 
    … …  
    14431473      ! 
    14441474   END SUBROUTINE updatee3t 
     1475#endif 
    14451476 
    14461477#else 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/NST/agrif_user.F90

    r13546 r14054  
    288288         CALL Agrif_Init_Variable(sshini_id, procname=agrif_initssh) 
    289289         CALL lbc_lnk( 'Agrif_Init_Domain', ssh(:,:,Kbb), 'T', 1. ) 
     290#if ! defined key_qco 
    290291         DO jk = 1, jpk 
    291292               e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb)  ) & 
    … …  
    293294                        &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    294295         END DO 
     296#endif 
    295297      ENDIF 
    296298 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DIA/diawri.F90

    r13497 r14054  
    1919   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output 
    2020   !!                 !                     change name of output variables in dia_wri_state 
     21   !!            4.0  ! 2020-10  (A. Nasser, S. Techene) add diagnostic for SWE 
    2122   !!---------------------------------------------------------------------- 
    2223 
    … …  
    4647   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
    4748   USE zdfmxl         ! mixed layer 
     49   USE zdfosm         ! mixed layer 
    4850   ! 
    4951   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    … …  
    118120      INTEGER ::   ji, jj, jk       ! dummy loop indices 
    119121      INTEGER ::   ikbot            ! local integer 
    120       REAL(wp)::   ze3 
    121122      REAL(wp)::   zztmp , zztmpx   ! local scalar 
    122123      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     124      REAL(wp)::   ze3 
    123125      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
    124126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
    … …  
    137139      CALL iom_put("e3u_0", e3u_0(:,:,:) ) 
    138140      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
     141      CALL iom_put("e3f_0", e3f_0(:,:,:) ) 
    139142      ! 
    140143      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t 
    … …  
    163166         CALL iom_put( "e3w" , z3d(:,:,:) ) 
    164167      ENDIF 
     168      IF ( iom_use("e3f") ) THEN                         ! time-varying e3f caution here at Kaa 
     169          DO jk = 1, jpk 
     170            z3d(:,:,jk) =  e3f(:,:,jk) 
     171         END DO 
     172         CALL iom_put( "e3f" , z3d(:,:,:) ) 
     173      ENDIF 
    165174 
    166175      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying) 
    167          CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) 
     176         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) ) 
    168177      ELSE 
    169178         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    170179      ENDIF 
    171180 
    172       IF( iom_use("wetdep") )   &                  ! wet depth 
    173          CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
     181      IF( iom_use("wetdep") )    CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )   ! wet depth 
     182          
     183#if defined key_qco 
     184      IF( iom_use("ht") )   CALL iom_put( "ht" , ht(:,:)     )   ! water column at t-point 
     185      IF( iom_use("hu") )   CALL iom_put( "hu" , hu(:,:,Kmm) )   ! water column at u-point 
     186      IF( iom_use("hv") )   CALL iom_put( "hv" , hv(:,:,Kmm) )   ! water column at v-point 
     187      IF( iom_use("hf") )   CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) )   ! water column at f-point (caution here at Naa) 
     188#endif 
    174189       
    175190      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
    … …  
    325340      ENDIF 
    326341      ! 
     342      IF ( iom_use("sKE") ) THEN                        ! surface kinetic energy at T point 
     343         z2d(:,:) = 0._wp 
     344         DO_2D( 0, 0, 0, 0 ) 
     345            z2d(ji,jj) = 0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  & 
     346               &                   + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm)  & 
     347               &                   + vv(ji,jj  ,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,1,Kmm)  &  
     348               &                   + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm)  )  & 
     349               &                 * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj) 
     350         END_2D 
     351         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
     352         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )    
     353      ENDIF 
     354      !     
     355      IF ( iom_use("sKEf") ) THEN                        ! surface kinetic energy at F point 
     356         z2d(:,:) = 0._wp                                ! CAUTION : only valid in SWE, not with bathymetry 
     357         DO_2D( 0, 0, 0, 0 ) 
     358            z2d(ji,jj) = 0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  & 
     359               &                   + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm)  & 
     360               &                   + vv(ji  ,jj,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji  ,jj) * e3v(ji  ,jj,1,Kmm)  &  
     361               &                   + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm)  )  & 
     362               &                 * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj) 
     363         END_2D 
     364         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     365         CALL iom_put( "sKEf", z2d )                      
     366      ENDIF 
     367      ! 
    327368      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
    328369 
    … …  
    424465       
    425466      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
     467       
     468      ! Output of vorticity terms 
     469      IF ( iom_use("relvor")    .OR. iom_use("plavor")    .OR.   & 
     470         & iom_use("relpotvor") .OR. iom_use("abspotvor") .OR.   & 
     471         & iom_use("Ens")                                        ) THEN 
     472         ! 
     473         z2d(:,:) = 0._wp  
     474         ze3 = 0._wp  
     475         DO_2D( 1, 0, 1, 0 ) 
     476            z2d(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    & 
     477            &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm)  ) * r1_e1e2f(ji,jj) 
     478         END_2D 
     479         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     480         CALL iom_put( "relvor", z2d )                  ! relative vorticity ( zeta )  
     481         ! 
     482         CALL iom_put( "plavor", ff_f )                 ! planetary vorticity ( f ) 
     483         ! 
     484         DO_2D( 1, 0, 1, 0 )   
     485            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     486              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     487            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     488            ELSE                      ;   ze3 = 0._wp 
     489            ENDIF 
     490            z2d(ji,jj) = ze3 * z2d(ji,jj)  
     491         END_2D 
     492         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     493         CALL iom_put( "relpotvor", z2d )                  ! relative potential vorticity (zeta/h) 
     494         ! 
     495         DO_2D( 1, 0, 1, 0 ) 
     496            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     497              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     498            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     499            ELSE                      ;   ze3 = 0._wp 
     500            ENDIF 
     501            z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)  
     502         END_2D 
     503         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     504         CALL iom_put( "abspotvor", z2d )                  ! absolute potential vorticity ( q ) 
     505         ! 
     506         DO_2D( 1, 0, 1, 0 )   
     507            z2d(ji,jj) = 0.5_wp * z2d(ji,jj)  * z2d(ji,jj)  
     508         END_2D 
     509         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     510         CALL iom_put( "Ens", z2d )                        ! potential enstrophy ( 1/2*q2 ) 
     511         ! 
     512      ENDIF 
    426513 
    427514      IF( ln_timing )   CALL timing_stop('dia_wri') 
    … …  
    9971084      !! 
    9981085      INTEGER :: inum, jk 
    999       REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept      ! 3D workspace !!st patch to use substitution 
     1086      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept       ! 3D workspace for qco substitution 
    10001087      !!---------------------------------------------------------------------- 
    10011088      !  
    … …  
    10761163         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
    10771164      ENDIF 
     1165      IF( ln_zdfosm ) THEN 
     1166         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)  )      ! now boundary-layer depth 
     1167         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)  )      ! now mixed-layer depth 
     1168         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask     )      ! w-level diffusion 
     1169         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask     )      ! now w-level viscosity 
     1170         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask     )      ! non-local t forcing 
     1171         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask     )      ! non-local s forcing 
     1172         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask     )      ! non-local u forcing 
     1173         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask     )      ! non-local v forcing 
     1174         IF( ln_osm_mle ) THEN 
     1175            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)  ) ! now transition-layer depth 
     1176         END IF 
     1177      ENDIF 
    10781178      ! 
    10791179      CALL iom_close( inum ) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/dom_oce.F90

    r14049 r14054  
    131131   ! 
    132132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2t , r1_e1e2t                !: associated metrics at t-point 
    133    REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2u , e2_e1u, r1_e1e2u        !: associated metrics at u-point 
    134    REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2v , e1_e2v, r1_e1e2v        !: associated metrics at v-point 
     133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2u , r1_e1e2u , e2_e1u       !: associated metrics at u-point 
     134   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2v , r1_e1e2v , e1_e2v       !: associated metrics at v-point 
    135135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point 
    136136   ! 
    … …  
    162162 
    163163   !                                                        !  reference depths of cells 
    164    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0  !: t- depth              [m] 
    165    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0  !: w- depth              [m] 
    166    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0  !: w- depth (sum of e3w) [m] 
     164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdept_0  !: t- depth              [m] 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdepw_0  !: w- depth              [m] 
     166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w_0  !: w- depth (sum of e3w) [m] 
    167167   !                                                        !  time-dependent depths of cells 
    168168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  gdept, gdepw   
    … …  
    205205 
    206206   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   ssmask, ssumask, ssvmask, ssfmask   !: surface mask at T-,U-, V- and F-pts 
    207    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask   !: land/ocean mask at T-, U-, V-, W- and F-pts 
    208    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    209  
     207   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   tmask, umask, vmask, wmask, fmask   !: land/ocean mask at T-, U-, V-, W- and F-pts 
     208   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   wumask, wvmask                      !: land/ocean mask at WU- and WV-pts 
     209#if defined key_qco    
     210   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   fe3mask                             !: land/ocean mask at F-pts for qco 
     211#endif 
    210212   !!---------------------------------------------------------------------- 
    211213   !! calendar variables 
    … …  
    306308         &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) ) 
    307309         ! 
    308 #if ! defined key_qco 
     310#if defined key_qco 
     311      ii = ii+1 
     312      ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,      & 
     313         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )              
     314#else 
    309315      ii = ii+1 
    310316      ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      & 
    … …  
    313319         ! 
    314320      ii = ii+1 
    315       ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  & 
    316          &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) ) 
    317          ! 
    318       ii = ii+1 
    319321      ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)    ,    hv_0(jpi,jpj)     , hf_0(jpi,jpj) ,       & 
    320322         &   r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) ,    r1_hv_0(jpi,jpj),   r1_hf_0(jpi,jpj) ,   STAT=ierr(ii)  ) 
    … …  
    323325      ii = ii+1 
    324326      ALLOCATE( ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
    325          &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
    326 #else 
    327       ii = ii+1 
    328       ALLOCATE(                    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
    329327         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
    330328#endif 
    … …  
    350348      ii = ii+1 
    351349      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) ) 
     350#if defined key_qco 
     351         ! 
     352      ii = ii+1 
     353      ALLOCATE( fe3mask(jpi,jpj,jpk) , STAT=ierr(ii) ) 
     354#endif 
    352355      ! 
    353356      dom_oce_alloc = MAXVAL(ierr) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/domain.F90

    r14049 r14054  
    1515   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default 
    1616   !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface 
    17    !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
     17   !!            4.1  !  2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    1818   !!---------------------------------------------------------------------- 
    1919    
    … …  
    2828   USE oce            ! ocean variables 
    2929   USE dom_oce        ! domain: ocean 
     30#if defined key_qco 
     31   USE domqco         ! quasi-eulerian 
     32#else 
     33   USE domvvl         ! variable volume 
     34#endif 
     35   USE sshwzv  , ONLY : ssh_init_rst   ! set initial ssh  
    3036   USE sbc_oce        ! surface boundary condition: ocean 
    3137   USE trc_oce        ! shared ocean & passive tracers variab 
    … …  
    3541   USE dommsk         ! domain: set the mask system 
    3642   USE domwri         ! domain: write the meshmask file 
    37 #if ! defined key_qco 
    38    USE domvvl         ! variable volume 
    39 #else 
    40    USE domqco          ! variable volume 
    41 #endif 
    4243   USE c1d            ! 1D configuration 
    4344   USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine) 
    44    USE wet_dry, ONLY : ll_wd 
    45    USE closea , ONLY : dom_clo ! closed seas 
     45   USE wet_dry , ONLY : ll_wd     ! wet & drying flag 
     46   USE closea  , ONLY : dom_clo   ! closed seas routine 
    4647   ! 
    4748   USE prtctl         ! Print control (prt_ctl_info routine) 
    … …  
    5051   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    5152   USE lib_mpp        ! distributed memory computing library 
     53   USE restart        ! only for lrst_oce 
    5254 
    5355   IMPLICIT NONE 
    … …  
    5860   PUBLIC   dom_tile     ! called by step.F90 
    5961 
     62   !! * Substitutions 
     63#  include "do_loop_substitute.h90" 
    6064   !!------------------------------------------------------------------------- 
    6165   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    8488      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices 
    8589      INTEGER ::   iconf = 0    ! local integers 
     90      REAL(wp)::   zrdt 
    8691      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"  
    8792      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    … …  
    121126         WRITE(numout,*)     '         cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg 
    122127      ENDIF 
    123       nn_wxios = 0 
    124       ln_xios_read = .FALSE. 
     128       
    125129      ! 
    126130      !           !==  Reference coordinate system  ==! 
    … …  
    143147      hv_0(:,:) = 0._wp 
    144148      hf_0(:,:) = 0._wp 
    145       DO jk = 1, jpk 
     149      DO jk = 1, jpkm1 
    146150         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    147151         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 
    148152         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 
    149          hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk) 
    150153      END DO 
     154      ! 
     155      DO jk = 1, jpkm1 
     156         hf_0(1:jpim1,:) = hf_0(1:jpim1,:) + e3f_0(1:jpim1,:,jk)*vmask(1:jpim1,:,jk)*vmask(2:jpi,:,jk) 
     157      END DO 
     158      CALL lbc_lnk('domain', hf_0, 'F', 1._wp) 
     159      ! 
     160      IF( lk_SWE ) THEN      ! SWE case redefine hf_0 
     161         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,1) * ssfmask(:,:) 
     162      ENDIF 
    151163      ! 
    152164      r1_ht_0(:,:) = ssmask (:,:) / ( ht_0(:,:) + 1._wp -  ssmask (:,:) ) 
    … …  
    154166      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp -  ssvmask(:,:) ) 
    155167      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp -  ssfmask(:,:) ) 
    156  
     168      ! 
     169      IF( ll_wd ) THEN       ! wet and drying (check ht_0 >= 0) 
     170         DO_2D( 1, 1, 1, 1 ) 
     171            IF( ht_0(ji,jj) < 0._wp .AND. ssmask(ji,jj) == 1._wp ) THEN 
     172               CALL ctl_stop( 'ssh_init_rst : ht_0 must be positive at potentially wet points' ) 
     173            ENDIF 
     174         END_2D 
     175      ENDIF 
     176      ! 
     177      !           !==  initialisation of time varying coordinate  ==! 
     178      ! 
     179      !                                 != ssh initialization 
     180      IF( .NOT.l_offline .AND. .NOT.l_SAS ) THEN 
     181         CALL ssh_init_rst( Kbb, Kmm, Kaa ) 
     182      ELSE 
     183         ssh(:,:,:) = 0._wp 
     184      ENDIF 
    157185      ! 
    158186#if defined key_qco 
    159       !           !==  initialisation of time varying coordinate  ==!  Quasi-Euerian coordinate case 
     187      !                                 != Quasi-Euerian coordinate case 
    160188      ! 
    161189      IF( .NOT.l_offline )   CALL dom_qco_init( Kbb, Kmm, Kaa ) 
    162       ! 
    163       IF( ln_linssh )        CALL ctl_stop('STOP','domain: key_qco and ln_linssh = T are incompatible') 
    164       ! 
    165190#else 
    166       !           !==  time varying part of coordinate system  ==! 
    167       ! 
    168       IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all 
     191      ! 
     192      IF( ln_linssh ) THEN              != Fix in time : set to the reference one for all 
    169193         ! 
    170194         DO jt = 1, jpt                         ! depth of t- and w-grid-points 
    … …  
    175199         ! 
    176200         DO jt = 1, jpt                         ! vertical scale factors 
    177             e3t(:,:,:,jt) =  e3t_0(:,:,:) 
    178             e3u(:,:,:,jt) =  e3u_0(:,:,:) 
    179             e3v(:,:,:,jt) =  e3v_0(:,:,:) 
    180             e3w(:,:,:,jt) =  e3w_0(:,:,:) 
     201            e3t (:,:,:,jt) =  e3t_0(:,:,:) 
     202            e3u (:,:,:,jt) =  e3u_0(:,:,:) 
     203            e3v (:,:,:,jt) =  e3v_0(:,:,:) 
     204            e3w (:,:,:,jt) =  e3w_0(:,:,:) 
    181205            e3uw(:,:,:,jt) = e3uw_0(:,:,:) 
    182206            e3vw(:,:,:,jt) = e3vw_0(:,:,:) 
    183207         END DO 
    184             e3f(:,:,:)    =  e3f_0(:,:,:) 
     208            e3f (:,:,:)    =  e3f_0(:,:,:) 
    185209         ! 
    186210         DO jt = 1, jpt                         ! water column thickness and its inverse 
    187             hu(:,:,jt)    =    hu_0(:,:) 
    188             hv(:,:,jt)    =    hv_0(:,:) 
     211               hu(:,:,jt) =    hu_0(:,:) 
     212               hv(:,:,jt) =    hv_0(:,:) 
    189213            r1_hu(:,:,jt) = r1_hu_0(:,:) 
    190214            r1_hv(:,:,jt) = r1_hv_0(:,:) 
    191215         END DO 
    192             ht(:,:) =    ht_0(:,:) 
    193          ! 
    194       ELSE                       != time varying : initialize before/now/after variables 
     216               ht   (:,:) =    ht_0(:,:) 
     217         ! 
     218      ELSE                              != Time varying : initialize before/now/after variables 
    195219         ! 
    196220         IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa ) 
    … …  
    373397      USE ioipsl 
    374398      !! 
    375       INTEGER  ::   ios   ! Local integer 
     399      INTEGER ::   ios   ! Local integer 
     400      REAL(wp)::   zrdt 
     401      !!---------------------------------------------------------------------- 
    376402      ! 
    377403      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
    … …  
    393419      ENDIF 
    394420      ! 
     421      !                       !=======================! 
     422      !                       !==  namelist namdom  ==! 
     423      !                       !=======================! 
     424      ! 
     425      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 
     426903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' ) 
     427      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
     428904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' ) 
     429      IF(lwm) WRITE( numond, namdom ) 
     430      ! 
     431#if defined key_agrif 
     432      IF( .NOT. Agrif_Root() ) THEN    ! AGRIF child, subdivide the Parent timestep 
     433         rn_Dt = Agrif_Parent (rn_Dt ) / Agrif_Rhot() 
     434      ENDIF 
     435#endif 
     436      ! 
     437      IF(lwp) THEN 
     438         WRITE(numout,*) 
     439         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain' 
     440         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
     441         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
     442         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt 
     443         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 
     444         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs 
     445      ENDIF 
     446      ! 
     447      ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3 
     448      rDt   = 2._wp * rn_Dt 
     449      r1_Dt = 1._wp / rDt 
     450      ! 
     451      IF( l_SAS .AND. .NOT.ln_linssh ) THEN 
     452         CALL ctl_warn( 'SAS requires linear ssh : force ln_linssh = T' ) 
     453         ln_linssh = .TRUE. 
     454      ENDIF 
     455      ! 
     456#if defined key_qco 
     457      IF( ln_linssh )   CALL ctl_stop( 'STOP','domain: key_qco and ln_linssh = T are incompatible' ) 
     458#endif 
     459      ! 
     460      !                       !=======================! 
     461      !                       !==  namelist namrun  ==! 
     462      !                       !=======================! 
    395463      ! 
    396464      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901) 
    … …  
    452520      nleapy = nn_leapy 
    453521      ninist = nn_istate 
     522      ! 
     523      !                                        !==  Set parameters for restart reading using xIOS  ==! 
     524      ! 
     525      IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
     526         lrxios = ln_xios_read .AND. ln_rstart 
     527         IF( nn_wxios > 0 )   lwxios = .TRUE.           !* set output file type for XIOS based on NEMO namelist 
     528         nxioso = nn_wxios 
     529      ENDIF 
     530      !                                        !==  Check consistency between ln_rstart and ln_1st_euler  ==!   (i.e. set l_1st_euler) 
    454531      l_1st_euler = ln_1st_euler 
    455       IF( .NOT. l_1st_euler .AND. .NOT. ln_rstart ) THEN 
     532      ! 
     533      IF( ln_rstart ) THEN                              !*  Restart case 
     534         ! 
     535         IF(lwp) WRITE(numout,*) 
     536         IF(lwp) WRITE(numout,*) '   open the restart file' 
     537         CALL rst_read_open                                              !- Open the restart file 
     538         ! 
     539         IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 ) THEN     !- Check time-step consistency and force Euler restart if changed 
     540            CALL iom_get( numror, 'rdt', zrdt ) 
     541            IF( zrdt /= rn_Dt ) THEN 
     542               IF(lwp) WRITE( numout,*) 
     543               IF(lwp) WRITE( numout,*) '   rn_Dt = ', rn_Dt,' not equal to the READ one rdt = ', zrdt 
     544               IF(lwp) WRITE( numout,*) 
     545               IF(lwp) WRITE( numout,*) '      ==>>>   forced euler first time-step' 
     546               l_1st_euler =  .TRUE. 
     547            ENDIF 
     548         ENDIF 
     549         ! 
     550         IF( .NOT.l_SAS .AND. iom_varid( numror, 'sshb', ldstop = .FALSE. ) <= 0 ) THEN   !- Check absence of one of the Kbb field (here sshb) 
     551            !                                                                             !  (any Kbb field is missing ==> all Kbb fields are missing)  
     552            IF( .NOT.l_1st_euler ) THEN 
     553               CALL ctl_warn('dom_nam : ssh at Kbb not found in restart files ',   & 
     554                  &                        'l_1st_euler forced to .true. and ' ,   & 
     555                  &                        'ssh(Kbb) = ssh(Kmm) '                  ) 
     556               l_1st_euler = .TRUE. 
     557            ENDIF 
     558         ENDIF 
     559      ELSEIF( .NOT.l_1st_euler ) THEN                   !*  Initialization case 
    456560         IF(lwp) WRITE(numout,*)   
    457561         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
    458562         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : l_1st_euler is forced to .true. '    
    459          l_1st_euler = .true. 
    460       ENDIF 
    461       !                             ! control of output frequency 
    462       IF( .NOT. ln_rst_list ) THEN     ! we use nn_stock 
     563         l_1st_euler = .TRUE. 
     564      ENDIF 
     565      ! 
     566      !                                        !==  control of output frequency  ==! 
     567      ! 
     568      IF( .NOT. ln_rst_list ) THEN   ! we use nn_stock 
    463569         IF( nn_stock == -1 )   CALL ctl_warn( 'nn_stock = -1 --> no restart will be done' ) 
    464570         IF( nn_stock == 0 .OR. nn_stock > nitend ) THEN 
    … …  
    479585      IF( Agrif_Root() ) THEN 
    480586         IF(lwp) WRITE(numout,*) 
    481          SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL 
     587         SELECT CASE ( nleapy )                !==  Choose calendar for IOIPSL  ==! 
    482588         CASE (  1 )  
    483589            CALL ioconf_calendar('gregorian') 
    … …  
    491597         END SELECT 
    492598      ENDIF 
    493  
    494       READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 
    495 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' ) 
    496       READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
    497 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' ) 
    498       IF(lwm) WRITE( numond, namdom ) 
    499       ! 
    500 #if defined key_agrif 
    501       IF( .NOT. Agrif_Root() ) THEN 
    502             rn_Dt = Agrif_Parent(rn_Dt) / Agrif_Rhot() 
    503       ENDIF 
    504 #endif 
    505       ! 
    506       IF(lwp) THEN 
    507          WRITE(numout,*) 
    508          WRITE(numout,*) '   Namelist : namdom   ---   space & time domain' 
    509          WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
    510          WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
    511          WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt 
    512          WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 
    513          WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs 
    514       ENDIF 
    515       ! 
    516       !! Initialise current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3 
    517       rDt  = 2._wp * rn_Dt 
    518       r1_Dt = 1._wp / rDt 
    519  
     599      ! 
     600      !                       !========================! 
     601      !                       !==  namelist namtile  ==! 
     602      !                       !========================! 
     603      ! 
    520604      READ  ( numnam_ref, namtile, IOSTAT = ios, ERR = 905 ) 
    521605905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtile in reference namelist' ) 
    … …  
    537621         ENDIF 
    538622      ENDIF 
    539  
    540       IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    541          lrxios = ln_xios_read.AND.ln_rstart 
    542 !set output file type for XIOS based on NEMO namelist  
    543          IF (nn_wxios > 0) lwxios = .TRUE.  
    544          nxioso = nn_wxios 
    545       ENDIF 
    546  
     623      ! 
    547624#if defined key_netcdf4 
    548       !                             ! NetCDF 4 case   ("key_netcdf4" defined) 
     625      !                       !=======================! 
     626      !                       !==  namelist namnc4  ==!   NetCDF 4 case   ("key_netcdf4" defined) 
     627      !                       !=======================! 
     628      ! 
    549629      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907) 
    550630907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namnc4 in reference namelist' ) 
    … …  
    555635      IF(lwp) THEN                        ! control print 
    556636         WRITE(numout,*) 
    557          WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters' 
     637         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters ("key_netcdf4" defined)' 
    558638         WRITE(numout,*) '      number of chunks in i-dimension             nn_nchunks_i = ', nn_nchunks_i 
    559639         WRITE(numout,*) '      number of chunks in j-dimension             nn_nchunks_j = ', nn_nchunks_j 
    … …  
    618698   SUBROUTINE domain_cfg( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 
    619699      !!---------------------------------------------------------------------- 
    620       !!                     ***  ROUTINE dom_nam  *** 
     700      !!                     ***  ROUTINE domain_cfg  *** 
    621701      !!                     
    622702      !! ** Purpose :   read the domain size in domain configuration file 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/dommsk.F90

    r13461 r14054  
    181181      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 
    182182      ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 
     183      IF( lk_SWE ) THEN      ! Shallow Water Eq. case : redefine ssfmask 
     184         DO_2D( 0,0 , 0,0 ) 
     185            ssfmask(ji,jj) = MAX(  ssmask(ji,jj+1), ssmask(ji+1,jj+1),  &  
     186               &                   ssmask(ji,jj  ), ssmask(ji+1,jj  )   ) 
     187         END_2D 
     188         CALL lbc_lnk( 'dommsk', ssfmask, 'F', 1.0_wp ) 
     189      ENDIF 
     190#if defined key_qco 
     191      fe3mask(:,:,:) = fmask(:,:,:) 
     192#endif 
    183193 
    184194      ! Interior domain mask  (used for global sum) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/domqco.F90

    r14049 r14054  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    10    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
    11    !!            4.x  !  2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate 
    12    !!---------------------------------------------------------------------- 
    13  
    14    !!---------------------------------------------------------------------- 
    15    !!   dom_qe_init   : define initial vertical scale factors, depths and column thickness 
    16    !!   dom_qe_r3c    : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points 
    17    !!       qe_rst_read : read/write restart file 
    18    !!   dom_qe_ctl    : Check the vvl options 
     10   !!            4.1  !  2019-08  (A. Coward, D. Storkey) add time level indices for prognostic variables 
     11   !!             -   !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*) 
     12   !!---------------------------------------------------------------------- 
     13 
     14   !!---------------------------------------------------------------------- 
     15   !!   dom_qco_init  : define initial vertical scale factors, depths and column thickness 
     16   !!   dom_qco_zgr   : Set ssh/h_0 ratio at t 
     17   !!   dom_qco_r3c   : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points 
     18   !!       qco_ctl   : Check the vvl options 
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce            ! ocean dynamics and tracers 
    … …  
    5555   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints 
    5656 
    57    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
    58  
    5957   !! * Substitutions 
    6058#  include "do_loop_substitute.h90" 
    … …  
    7977      !! 
    8078      !!---------------------------------------------------------------------- 
    81       INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     79      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices 
     80      !!---------------------------------------------------------------------- 
    8281      ! 
    8382      IF(lwp) WRITE(numout,*) 
    … …  
    8584      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    8685      ! 
    87       CALL dom_qco_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    88       ! 
    89       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    90       CALL qe_rst_read( nit000, Kbb, Kmm ) 
    91       ! 
    92       CALL dom_qco_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
     86      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer) 
     87      ! 
     88      CALL dom_qco_zgr( Kbb, Kmm )            ! interpolation scale factor, depth and water column 
     89      ! 
     90#if defined key_agrif 
     91      ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a 
     92      ! problem for the restartability...) 
     93      r3t(:,:,Kaa) = r3t(:,:,Kmm) 
     94      r3u(:,:,Kaa) = r3u(:,:,Kmm) 
     95      r3v(:,:,Kaa) = r3v(:,:,Kmm) 
     96#endif 
    9397      ! 
    9498   END SUBROUTINE dom_qco_init 
    9599 
    96100 
    97    SUBROUTINE dom_qco_zgr(Kbb, Kmm, Kaa) 
     101   SUBROUTINE dom_qco_zgr( Kbb, Kmm ) 
    98102      !!---------------------------------------------------------------------- 
    99103      !!                ***  ROUTINE dom_qco_init  *** 
    100104      !! 
    101       !! ** Purpose :  Initialization of all ssh. to h._0 ratio 
    102       !! 
    103       !! ** Method  :  - interpolate scale factors 
    104       !! 
    105       !! ** Action  : - r3(t/u/v)_b 
    106       !!              - r3(t/u/v/f)_n 
    107       !! 
    108       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    109       !!---------------------------------------------------------------------- 
    110       INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     105      !! ** Purpose :  Initialization of all r3. = ssh./h._0 ratios 
     106      !! 
     107      !! ** Method  :  Call domqco using Kbb and Kmm 
     108      !!               NB: dom_qco_zgr is called by dom_qco_init it uses ssh from ssh_init  
     109      !! 
     110      !! ** Action  : - r3(t/u/v)(Kbb) 
     111      !!              - r3(t/u/v/f)(Kmm) 
     112      !!---------------------------------------------------------------------- 
     113      INTEGER, INTENT(in) ::   Kbb, Kmm   ! time level indices 
    111114      !!---------------------------------------------------------------------- 
    112115      ! 
    113116      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    114117      !                                ! Horizontal interpolation of e3t 
    115       CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) ) 
     118      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           ) 
    116119      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 
    117120      ! 
    … …  
    143146      !                                      !==  ratio at u-,v-point  ==! 
    144147      ! 
    145       IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging) 
     148!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging) 
     149#if ! defined key_qcoTest_FluxForm 
     150      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
    146151         DO_2D( 0, 0, 0, 0 ) 
    147152            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  & 
    … …  
    150155               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 
    151156         END_2D 
    152       ELSE                                         !- Flux Form   (simple averaging) 
     157!!st      ELSE                                         !- Flux Form   (simple averaging) 
     158#else 
    153159         DO_2D( 0, 0, 0, 0 ) 
    154             pr3u(ji,jj) = 0.5_wp * (  pssh(ji  ,jj) + pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) 
    155             pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj  ) + pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) 
     160            pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj) 
     161            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj) 
    156162         END_2D 
    157       ENDIF 
     163!!st      ENDIF 
     164#endif          
    158165      ! 
    159166      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only 
    … …  
    163170      ELSE                                   !==  ratio at f-point  ==! 
    164171         ! 
    165          IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging) 
    166             DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
     172!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging) 
     173#if ! defined key_qcoTest_FluxForm 
     174         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
     175 
     176            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
    167177               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  & 
    168178                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  & 
    … …  
    170180                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 
    171181            END_2D 
    172          ELSE                                      !- Flux Form   (simple averaging) 
    173             DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
    174                pr3f(ji,jj) = 0.25_wp * (  pssh(ji  ,jj  ) + pssh(ji+1,jj  )  & 
    175                   &                     + pssh(ji  ,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) 
     182!!st         ELSE                                      !- Flux Form   (simple averaging) 
     183#else 
     184            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
     185               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  & 
     186                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) 
    176187            END_2D 
    177          ENDIF 
     188!!st         ENDIF 
     189#endif 
    178190         !                                                 ! lbc on ratio at u-,v-,f-points 
    179191         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp ) 
    … …  
    184196 
    185197 
    186    SUBROUTINE qe_rst_read( kt, Kbb, Kmm ) 
     198   SUBROUTINE qco_ctl 
    187199      !!--------------------------------------------------------------------- 
    188       !!                   ***  ROUTINE qe_rst_read  *** 
    189       !! 
    190       !! ** Purpose :   Read ssh in restart file 
    191       !! 
    192       !! ** Method  :   use of IOM library 
    193       !!                if the restart does not contain ssh, 
    194       !!                it is set to the _0 values. 
    195       !!---------------------------------------------------------------------- 
    196       INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
    197       INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices 
    198       ! 
    199       INTEGER ::   ji, jj, jk 
    200       INTEGER ::   id1, id2     ! local integers 
    201       !!---------------------------------------------------------------------- 
    202       ! 
    203          IF( ln_rstart ) THEN                   !* Read the restart file 
    204             CALL rst_read_open                  !  open the restart file if necessary 
    205             ! 
    206             id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 
    207             id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 
    208             ! 
    209             !                             ! --------- ! 
    210             !                             ! all cases ! 
    211             !                             ! --------- ! 
    212             ! 
    213             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    214                CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb)    ) 
    215                CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm)    ) 
    216                ! needed to restart if land processor not computed 
    217                IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 
    218                WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh when it was required on e3 
    219                   ssh(:,:,Kmm) = 0._wp 
    220                   ssh(:,:,Kbb) = 0._wp 
    221                END WHERE 
    222                IF( l_1st_euler ) THEN 
    223                   ssh(:,:,Kbb) = ssh(:,:,Kmm) 
    224                ENDIF 
    225             ELSE IF( id1 > 0 ) THEN 
    226                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 
    227                IF(lwp) write(numout,*) 'sshn set equal to sshb.' 
    228                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    229                CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) ) 
    230                ssh(:,:,Kmm) = ssh(:,:,Kbb) 
    231                l_1st_euler = .TRUE. 
    232             ELSE IF( id2 > 0 ) THEN 
    233                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 
    234                IF(lwp) write(numout,*) 'sshb set equal to sshn.' 
    235                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    236                CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm) ) 
    237                ssh(:,:,Kbb) = ssh(:,:,Kmm) 
    238                l_1st_euler = .TRUE. 
    239             ELSE 
    240                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 
    241                IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 
    242                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    243                ssh(:,:,:) = 0._wp 
    244                l_1st_euler = .TRUE. 
    245             ENDIF 
    246             ! 
    247          ELSE                                   !* Initialize at "rest" 
    248             ! 
    249             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
    250                ! 
    251                IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case 
    252                   CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    253                   ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    254                   ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb) 
    255                   uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
    256                   vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    257                ELSE                                  ! if not test case 
    258                   ssh(:,:,Kmm) = -ssh_ref 
    259                   ssh(:,:,Kbb) = -ssh_ref 
    260                   ! 
    261                   DO_2D( 1, 1, 1, 1 ) 
    262                      IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    263                         ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    264                         ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    265                      ENDIF 
    266                   END_2D 
    267                ENDIF 
    268  
    269                DO ji = 1, jpi 
    270                   DO jj = 1, jpj 
    271                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    272                        CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' ) 
    273                      ENDIF 
    274                   END DO 
    275                END DO 
    276                ! 
    277             ELSE 
    278                ! 
    279                ! Just to read set ssh in fact, called latter once vertical grid 
    280                ! is set up: 
    281 !               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    282 !               ! 
    283                 ssh(:,:,:) = 0._wp 
    284                ! 
    285             ENDIF           ! end of ll_wd edits 
    286             ! 
    287          ENDIF 
    288       ! 
    289    END SUBROUTINE qe_rst_read 
    290  
    291  
    292    SUBROUTINE dom_qco_ctl 
    293       !!--------------------------------------------------------------------- 
    294       !!                  ***  ROUTINE dom_qco_ctl  *** 
     200      !!                  ***  ROUTINE qco_ctl  *** 
    295201      !! 
    296202      !! ** Purpose :   Control the consistency between namelist options 
    … …  
    312218      IF(lwp) THEN                    ! Namelist print 
    313219         WRITE(numout,*) 
    314          WRITE(numout,*) 'dom_qco_ctl : choice/control of the variable vertical coordinate' 
    315          WRITE(numout,*) '~~~~~~~~~~~' 
     220         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate' 
     221         WRITE(numout,*) '~~~~~~~~' 
    316222         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate' 
    317223         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar 
    … …  
    357263#endif 
    358264      ! 
    359    END SUBROUTINE dom_qco_ctl 
     265   END SUBROUTINE qco_ctl 
    360266 
    361267   !!====================================================================== 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/domvvl.F90

    r14049 r14054  
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1010   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
    11    !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
     11   !!             -   ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    1212   !!---------------------------------------------------------------------- 
    1313 
    … …  
    766766      !! ** Purpose :   Read or write VVL file in restart file 
    767767      !! 
    768       !! ** Method  :   use of IOM library 
    769       !!                if the restart does not contain vertical scale factors, 
    770       !!                they are set to the _0 values 
    771       !!                if the restart does not contain vertical scale factors increments (z_tilde), 
    772       !!                they are set to 0. 
     768      !! ** Method  : * restart comes from a linear ssh simulation : 
     769      !!                   an attempt to read e3t_n stops simulation 
     770      !!              * restart comes from a z-star, z-tilde, or layer : 
     771      !!                   read e3t_n and e3t_b 
     772      !!              * restart comes from a z-star : 
     773      !!                   set tilde_e3t_n, tilde_e3t_n, and hdiv_lf to 0 
     774      !!              * restart comes from layer : 
     775      !!                   read tilde_e3t_n and tilde_e3t_b 
     776      !!                   set hdiv_lf to 0 
     777      !!              * restart comes from a z-tilde: 
     778      !!                   read tilde_e3t_n, tilde_e3t_b, and hdiv_lf 
     779      !! 
     780      !!              NB: if l_1st_euler = T (ln_1st_euler or ssh_b not found) 
     781      !!                   Kbb fields set to Kmm ones 
    773782      !!---------------------------------------------------------------------- 
    774783      INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
    … …  
    776785      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag 
    777786      ! 
    778       INTEGER ::   ji, jj, jk 
    779       INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    780       !!---------------------------------------------------------------------- 
    781       ! 
    782       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    783          !                                   ! =============== 
    784          IF( ln_rstart ) THEN                   !* Read the restart file 
    785             CALL rst_read_open                  !  open the restart file if necessary 
    786             CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm)    ) 
     787      INTEGER ::   ji, jj, jk      ! dummy loop indices 
     788      INTEGER ::   id3, id4, id5   ! local integers 
     789      !!---------------------------------------------------------------------- 
     790      ! 
     791      !                                      !=====================! 
     792      IF( TRIM(cdrw) == 'READ' ) THEN        !  Read / initialise  ! 
     793         !                                   !=====================! 
     794         ! 
     795         IF( ln_rstart ) THEN                   !==  Read the restart file  ==! 
    787796            ! 
    788             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
    789             id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 
    790             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
     797            CALL rst_read_open                                          !*  open the restart file if necessary 
     798            !                                         ! --------- ! 
     799            !                                         ! all cases ! 
     800            !                                         ! --------- ! 
     801            ! 
     802            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )  !*  check presence 
    791803            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    792             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     804            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. ) 
    793805            ! 
    794             !                             ! --------- ! 
    795             !                             ! all cases ! 
    796             !                             ! --------- ! 
    797             ! 
    798             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     806            !                                                           !*  scale factors 
     807            IF(lwp) WRITE(numout,*)    '          Kmm scale factor read in the restart file' 
     808            CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) ) 
     809            WHERE ( tmask(:,:,:) == 0.0_wp )  
     810               e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
     811            END WHERE 
     812            IF( l_1st_euler ) THEN                       ! euler 
     813               IF(lwp) WRITE(numout,*) '          Euler first time step : e3t(Kbb) = e3t(Kmm)' 
     814               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     815            ELSE                                         ! leap frog 
     816               IF(lwp) WRITE(numout,*) '          Kbb scale factor read in the restart file' 
    799817               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) ) 
    800                CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) ) 
    801                ! needed to restart if land processor not computed  
    802                IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    803818               WHERE ( tmask(:,:,:) == 0.0_wp )  
    804                   e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
    805819                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    806820               END WHERE 
    807                IF( l_1st_euler ) THEN 
    808                   e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    809                ENDIF 
    810             ELSE IF( id1 > 0 ) THEN 
    811                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    812                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    813                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    814                CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) ) 
    815                e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    816                l_1st_euler = .true. 
    817             ELSE IF( id2 > 0 ) THEN 
    818                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    819                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    820                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    821                CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) ) 
    822                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    823                l_1st_euler = .true. 
    824             ELSE 
    825                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    826                IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    827                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    828                DO jk = 1, jpk 
    829                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
    830                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    831                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    832                END DO 
    833                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    834                l_1st_euler = .true. 
    835821            ENDIF 
    836             !                             ! ----------- ! 
    837             IF( ln_vvl_zstar ) THEN       ! z_star case ! 
    838                !                          ! ----------- ! 
     822            !                                         ! ------------ ! 
     823            IF( ln_vvl_zstar ) THEN                   ! z_star case ! 
     824               !                                      ! ------------ ! 
    839825               IF( MIN( id3, id4 ) > 0 ) THEN 
    840826                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 
    841827               ENDIF 
    842                !                          ! ----------------------- ! 
    843             ELSE                          ! z_tilde and layer cases ! 
    844                !                          ! ----------------------- ! 
    845                IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
    846                   CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
     828               !                                      ! ------------------------ ! 
     829            ELSE                                      !  z_tilde and layer cases ! 
     830               !                                      ! ------------------------ ! 
     831               ! 
     832               IF( id4 > 0 ) THEN                                       !*  scale factor increments 
     833                  IF(lwp) WRITE(numout,*)    '          Kmm scale factor increments read in the restart file' 
    847834                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 
    848                ELSE                            ! one at least array is missing 
     835                  IF( l_1st_euler ) THEN                 ! euler 
     836                     IF(lwp) WRITE(numout,*) '          Euler first time step : tilde_e3t(Kbb) = tilde_e3t(Kmm)' 
     837                     tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 
     838                  ELSE                                   ! leap frog 
     839                     IF(lwp) WRITE(numout,*) '          Kbb scale factor increments read in the restart file' 
     840                     CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
     841                  ENDIF 
     842               ELSE  
    849843                  tilde_e3t_b(:,:,:) = 0.0_wp 
    850844                  tilde_e3t_n(:,:,:) = 0.0_wp 
    851845               ENDIF 
    852                !                          ! ------------ ! 
    853                IF( ln_vvl_ztilde ) THEN   ! z_tilde case ! 
    854                   !                       ! ------------ ! 
     846               !                                      ! ------------ ! 
     847               IF( ln_vvl_ztilde ) THEN               ! z_tilde case ! 
     848                  !                                   ! ------------ ! 
    855849                  IF( id5 > 0 ) THEN  ! required array exists 
    856850                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    857851                  ELSE                ! array is missing 
    858                      hdiv_lf(:,:,:) = 0.0_wp 
     852                     hdiv_lf(:,:,:) = 0.0_wp  
    859853                  ENDIF 
    860854               ENDIF 
    861855            ENDIF 
    862856            ! 
    863          ELSE                                   !* Initialize at "rest" 
     857         ELSE                                   !==  Initialize at "rest" with ssh  ==! 
    864858            ! 
    865  
    866             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential  
    867                ! 
    868                IF( cn_cfg == 'wad' ) THEN 
    869                   ! Wetting and drying test case 
    870                   CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    871                   ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    872                   ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
    873                   uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    874                   vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    875                ELSE 
    876                   ! if not test case 
    877                   ssh(:,:,Kmm) = -ssh_ref 
    878                   ssh(:,:,Kbb) = -ssh_ref 
    879  
    880                   DO_2D( 1, 1, 1, 1 ) 
    881                      IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    882                         ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    883                         ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    884                      ENDIF 
    885                   END_2D 
    886                ENDIF !If test case else 
    887  
    888                ! Adjust vertical metrics for all wad 
    889                DO jk = 1, jpk 
    890                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) & 
    891                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    892                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    893                END DO 
    894                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    895  
    896                DO_2D( 1, 1, 1, 1 ) 
    897                   IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    898                      CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    899                   ENDIF 
    900                END_2D 
    901                ! 
    902             ELSE 
    903                ! 
    904                ! Just to read set ssh in fact, called latter once vertical grid 
    905                ! is set up: 
    906 !               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    907 !               ! 
    908 !               DO jk=1,jpk 
    909 !                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    910 !                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 
    911 !               END DO 
    912 !               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    913                 ssh(:,:,Kmm)=0._wp 
    914                 e3t(:,:,:,Kmm)=e3t_0(:,:,:) 
    915                 e3t(:,:,:,Kbb)=e3t_0(:,:,:) 
    916                ! 
    917             END IF           ! end of ll_wd edits 
    918  
     859            DO jk = 1, jpk 
     860               e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) ) 
     861            END DO 
     862            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     863            ! 
    919864            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    920865               tilde_e3t_b(:,:,:) = 0._wp 
    921866               tilde_e3t_n(:,:,:) = 0._wp 
    922867               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 
    923             END IF 
     868            ENDIF 
    924869         ENDIF 
    925          ! 
    926       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    927          !                                   ! =================== 
     870         !                                       !=======================! 
     871      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN       !  Create restart file  ! 
     872         !                                       !=======================! 
     873         ! 
    928874         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 
    929875         !                                           ! --------- ! 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/domzgr_substitute.h90

    r13237 r14054  
    1515#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1._wp+r3u(i,j,t)*umask(i,j,k))) 
    1616#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k))) 
    17 #   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k))) 
     17#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fe3mask(i,j,k))) 
     18#   define  e3f_vor(i,j,k) (e3f_0vor(i,j,k)*(1._wp+r3f(i,j)*fe3mask(i,j,k))) 
    1819#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t))) 
    1920#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t))) 
    2021#   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1._wp+r3v(i,j,t))) 
    21 #   define  ht(i,j)        (ht_0(i,j)+ssh(i,j,Kmm)) 
     22#   define  ht(i,j)        (ht_0(i,j)*(1._wp+r3t(i,j,Kmm))) 
    2223#   define  hu(i,j,t)      (hu_0(i,j)*(1._wp+r3u(i,j,t))) 
    2324#   define  hv(i,j,t)      (hv_0(i,j)*(1._wp+r3v(i,j,t))) 
    … …  
    2930#endif 
    3031!!---------------------------------------------------------------------- 
     32!!#   define  e3t_f(i,j,k)   (e3t_0(i,j,k)*(1._wp+r3t_f(i,j)*tmask(i,j,k))) 
     33!!#   define  e3u_f(i,j,k)   (e3u_0(i,j,k)*(1._wp+r3u_f(i,j)*umask(i,j,k))) 
     34!!#   define  e3v_f(i,j,k)   (e3v_0(i,j,k)*(1._wp+r3v_f(i,j)*vmask(i,j,k))) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/istate.F90

    r13295 r14054  
    4242   PRIVATE 
    4343 
    44    PUBLIC   istate_init   ! routine called by step.F90 
     44   PUBLIC   istate_init   ! routine called by nemogcm.F90 
    4545 
    4646   !! * Substitutions 
    … …  
    5959      !!  
    6060      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
     61      !! 
     62      !! ** Method  :    
    6163      !!---------------------------------------------------------------------- 
    6264      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices 
    6365      ! 
    6466      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    65       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute 
     67      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table for qco substitute 
    6668!!gm see comment further down 
    6769      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace 
    … …  
    7375      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    7476 
    75 !!gm  Why not include in the first call of dta_tsd ?   
    76 !!gm  probably associated with the use of internal damping... 
    7777       CALL dta_tsd_init        ! Initialisation of T & S input data 
    78 !!gm to be moved in usrdef of C1D case 
     78 
    7979!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data 
    80 !!gm 
    8180 
    82       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    83       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    84       ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    85       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
     81      rhd  (:,:,:      ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
     82      rn2b (:,:,:      ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
     83      ts   (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
     84      rab_b(:,:,:,:    ) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    8685#if defined key_agrif 
    8786      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization 
    … …  
    9695         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent 
    9796         ! 
    98          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)  
    99          ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
    100          uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    101          vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     97         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
     98!!st 
     99!!st need for a recent agrif version to be displaced toward ssh_init_rst with agrif_istate_ssh 
     100         ssh(:,:,    Kmm) = ssh(:,:    ,Kbb) 
     101!!st end 
     102         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     103         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    102104      ELSE 
    103105#endif 
    … …  
    117119            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000 
    118120            ! 
    119             ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest 
    120             uu  (:,:,:,Kbb) = 0._wp 
    121             vv  (:,:,:,Kbb) = 0._wp   
     121            uu (:,:,:,Kbb) = 0._wp 
     122            vv (:,:,:,Kbb) = 0._wp   
    122123            ! 
    123             IF( ll_wd ) THEN 
    124                ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
    125                ! 
    126                ! Apply minimum wetdepth criterion 
    127                ! 
    128                DO_2D( 1, 1, 1, 1 ) 
    129                   IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
    130                      ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    131                   ENDIF 
    132                END_2D 
    133             ENDIF  
    134              ! 
    135124         ELSE                                 ! user defined initial T and S 
    136125            DO jk = 1, jpk 
    137126               zgdept(:,:,jk) = gdept(:,:,jk,Kbb) 
    138127            END DO 
    139             CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
     128            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) ) 
    140129         ENDIF 
    141          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    142          ssh (:,:,Kmm)     = ssh(:,:,Kbb)    
    143          uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    144          vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    145  
    146 !!gm POTENTIAL BUG : 
    147 !!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    148 !!             as well as gdept_ and gdepw_....   !!!!!  
    149 !!      ===>>>>   probably a call to domvvl initialisation here.... 
    150  
     130         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     131         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     132         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    151133 
    152134         ! 
    153 !!gm to be moved in usrdef of C1D case 
    154 !         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 
    155 !            ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 
    156 !            CALL dta_uvd( nit000, zuvd ) 
    157 !            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 
    158 !            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 
    159 !            DEALLOCATE( zuvd ) 
    160 !         ENDIF 
     135!!gm ==>>>  to be moved in usrdef_istate of C1D case  
     136         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 
     137            ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 
     138            CALL dta_uvd( nit000, Kbb, zuvd ) 
     139            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 
     140            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 
     141            DEALLOCATE( zuvd ) 
     142         ENDIF 
    161143         ! 
    162 !!gm This is to be changed !!!! 
    163 !         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here 
    164 !         IF( .NOT.ln_linssh ) THEN 
    165 !            DO jk = 1, jpk 
    166 !               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    167 !            END DO 
    168 !         ENDIF 
    169 !!gm  
    170144         !  
    171145      ENDIF  

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DOM/phycst.F90

    r12489 r14054  
    6666   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos 
    6767   REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi 
     68    
    6869   !!---------------------------------------------------------------------- 
    6970   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/dynadv.F90

    r12377 r14054  
    127127      IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' ) 
    128128      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 
    129  
     129#if defined key_qcoTest_FluxForm 
     130      IF( ln_dynadv_vec  ) THEN CALL ctl_stop( 'STOP', 'key_qcoTest_FluxForm requires flux form advection' ) 
     131#endif 
    130132 
    131133      IF(lwp) THEN                    ! Print the choice 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/dynatf_qco.F90

    r13295 r14054  
    1 MODULE dynatfqco 
     1MODULE dynatf_qco 
    22   !!========================================================================= 
    3    !!                       ***  MODULE  dynatfqco  *** 
     3   !!                       ***  MODULE  dynatf_qco  *** 
    44   !! Ocean dynamics: time filtering 
    55   !!========================================================================= 
    … …  
    5050   USE prtctl         ! Print control 
    5151   USE timing         ! Timing 
    52 #if defined key_agrif 
    53    USE agrif_oce_interp 
    54 #endif 
    5552 
    5653   IMPLICIT NONE 
    … …  
    199196      ! JC: Would be more clever to swap variables than to make a full vertical 
    200197      ! integration 
    201       ! 
     198      ! CAUTION : calculation need to be done in the same way than see GM   
    202199      uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
    203       uu_b(:,:,Kmm) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
     200      uu_b(:,:,Kmm) = (e3u_0(:,:,1) * ( 1._wp + r3u_f(:,:) * umask(:,:,1) )) * puu(:,:,1,Kmm) * umask(:,:,1) 
    204201      vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
    205       vv_b(:,:,Kmm) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
     202      vv_b(:,:,Kmm) = (e3v_0(:,:,1) * ( 1._wp + r3v_f(:,:) * vmask(:,:,1))) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    206203      DO jk = 2, jpkm1 
    207204         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
    208          uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
     205         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + (e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
    209206         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
    210          vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
     207         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + (e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    211208      END DO 
    212209      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa) 
    213210      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa) 
    214       uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
    215       vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 
     211      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * (r1_hu_0(:,:)/( 1._wp + r3u_f(:,:) )) 
     212      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * (r1_hv_0(:,:)/( 1._wp + r3v_f(:,:) )) 
    216213      ! 
    217214      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents 
    … …  
    235232 
    236233   !!========================================================================= 
    237 END MODULE dynatfqco 
     234END MODULE dynatf_qco 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/dynldf_lap_blp.F90

    r13497 r14054  
    55   !!====================================================================== 
    66   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
     7   !!           4.0  ! 2020-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor  
    78   !!---------------------------------------------------------------------- 
    89 
    … …  
    1920   USE in_out_manager ! I/O manager 
    2021   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    21  
     22   USE lib_mpp 
     23    
    2224   IMPLICIT NONE 
    2325   PRIVATE 
    … …  
    4749      !! 
    4850      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 
     51      !! 
     52      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/  
    4953      !!---------------------------------------------------------------------- 
    5054      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    … …  
    5761      REAL(wp) ::   zsign        ! local scalars 
    5862      REAL(wp) ::   zua, zva     ! local scalars 
    59       REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
     63      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zcur, zdiv 
     64      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms 
    6065      !!---------------------------------------------------------------------- 
    6166      ! 
    … …  
    7075      ENDIF 
    7176      ! 
    72       !                                                ! =============== 
    73       DO jk = 1, jpkm1                                 ! Horizontal slab 
    74          !                                             ! =============== 
    75          DO_2D( 0, 1, 0, 1 ) 
    76             !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    77             zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask 
    78                &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
    79                &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    80             !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    81             zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask 
    82                &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
    83                &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
    84          END_2D 
     77      SELECT CASE( nn_dynldf_typ )   
     78      !               
     79      CASE ( np_typ_rot )       !==  Vorticity-Divergence operator  ==! 
    8580         ! 
    86          DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div ) 
    87             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use 
    88                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
    89                &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      ) 
     81         ALLOCATE( zcur(jpi,jpj) , zdiv(jpi,jpj) ) 
     82         ! 
     83         DO jk = 1, jpkm1                                 ! Horizontal slab 
     84            ! 
     85            DO_2D( 0, 1, 0, 1 ) 
     86               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
     87               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask 
     88                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     89                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
     90               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
     91               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask 
     92                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
     93                  &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
     94            END_2D 
     95            ! 
     96            DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div ) 
     97               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use 
     98                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
     99                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      ) 
    90100               ! 
    91             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use 
    92                &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
    93                &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      ) 
    94          END_2D 
    95          !                                             ! =============== 
    96       END DO                                           !   End of slab 
    97       !                                                ! =============== 
     101               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use 
     102                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
     103                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      ) 
     104            END_2D 
     105            ! 
     106         END DO                                           !   End of slab 
     107         ! 
     108         DEALLOCATE( zcur , zdiv ) 
     109         ! 
     110      CASE ( np_typ_sym )       !==  Symmetric operator  ==! 
     111         ! 
     112         ALLOCATE( zten(jpi,jpj) , zshe(jpi,jpj) ) 
     113         ! 
     114         DO jk = 1, jpkm1                                 ! Horizontal slab 
     115            ! 
     116            DO_2D( 0, 1, 0, 1 ) 
     117               !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask 
     118               zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              & 
     119                  &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             & 
     120                  &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   & 
     121                  &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             & 
     122                  &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   )  
     123               !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask 
     124               zten(ji,jj)    = ahmt(ji,jj,jk)                                                       & 
     125                  &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         & 
     126                  &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   & 
     127                  &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         & 
     128                  &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )    
     129            END_2D 
     130            ! 
     131            DO_2D( 0, 0, 0, 0 ) 
     132               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               & 
     133                  &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       & 
     134                  &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                     
     135                  &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     136                  &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )    
     137               ! 
     138               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               & 
     139                  &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     140                  &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     & 
     141                  &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       & 
     142                  &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) ) 
     143               ! 
     144            END_2D 
     145            ! 
     146         END DO 
     147         ! 
     148         DEALLOCATE( zten , zshe ) 
     149         ! 
     150      END SELECT 
    98151      ! 
    99152   END SUBROUTINE dyn_ldf_lap 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/dynspg_ts.F90

    r14049 r14054  
    306306      ENDIF 
    307307      ! 
    308       !                                   !=  Add atmospheric pressure forcing  =! 
    309       !                                   !  ----------------------------------  ! 
    310       IF( ln_bt_fw ) THEN                        ! Add wind forcing 
     308      !                                   !=  Add wind forcing  =! 
     309      !                                   !  ------------------  ! 
     310      IF( ln_bt_fw ) THEN 
    311311         DO_2D( 0, 0, 0, 0 ) 
    312312            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 
    … …  
    386386      ! 
    387387      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 
    388          zhup2_e(:,:) = hu(:,:,Kmm) 
    389          zhvp2_e(:,:) = hv(:,:,Kmm) 
    390          zhtp2_e(:,:) = ht(:,:) 
    391       ENDIF 
    392       ! 
    393       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    394          sshn_e(:,:) =    pssh(:,:,Kmm)             
     388         zhup2_e(:,:) = hu_0(:,:) 
     389         zhvp2_e(:,:) = hv_0(:,:) 
     390         zhtp2_e(:,:) = ht_0(:,:) 
     391      ENDIF 
     392      ! 
     393      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields                     
     394         sshn_e(:,:) =    pssh (:,:,Kmm)             
    395395         un_e  (:,:) =    puu_b(:,:,Kmm)             
    396396         vn_e  (:,:) =    pvv_b(:,:,Kmm) 
    … …  
    401401         hvr_e (:,:) = r1_hv(:,:,Kmm) 
    402402      ELSE                                ! CENTRED integration: start from BEFORE fields 
    403          sshn_e(:,:) =    pssh(:,:,Kbb) 
     403         sshn_e(:,:) =    pssh (:,:,Kbb) 
    404404         un_e  (:,:) =    puu_b(:,:,Kbb)          
    405405         vn_e  (:,:) =    pvv_b(:,:,Kbb) 
    … …  
    412412      ! 
    413413      ! Initialize sums: 
    414       puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)           
    415       pvv_b  (:,:,Kaa) = 0._wp 
     414      puu_b (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     415      pvv_b (:,:,Kaa) = 0._wp 
    416416      pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level 
    417       un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    418       vn_adv(:,:) = 0._wp 
     417      un_adv(:,:)     = 0._wp       ! Sum for now transport issued from ts loop 
     418      vn_adv(:,:)     = 0._wp 
    419419      ! 
    420420      IF( ln_wd_dl ) THEN 
    … …  
    475475            ! 
    476476            !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk) 
     477#if defined key_qcoTest_FluxForm 
     478            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     479            DO_2D( 1, 1, 1, 0 )   ! not jpi-column 
     480               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     481            END_2D 
     482            DO_2D( 1, 0, 1, 1 ) 
     483               zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     484            END_2D 
     485#else 
     486            !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
    477487            DO_2D( 1, 1, 1, 0 )   ! not jpi-column 
    478488               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        & 
    … …  
    485495                    &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
    486496            END_2D 
     497#endif                
    487498            ! 
    488499         ENDIF 
    … …  
    540551         !   
    541552         ! Sea Surface Height at u-,v-points (vvl case only) 
    542          IF( .NOT.ln_linssh ) THEN                                 
     553         IF( .NOT.ln_linssh ) THEN 
     554#if defined key_qcoTest_FluxForm 
     555            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     556            DO_2D( 1, 1, 1, 0 ) 
     557               zsshu_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     558            END_2D 
     559            DO_2D( 1, 0, 1, 1 ) 
     560               zsshv_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     561            END_2D 
     562#else 
    543563            DO_2D( 0, 0, 0, 0 ) 
    544                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    & 
    545                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
    546                   &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) ) 
    547                zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    & 
    548                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
    549                   &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) ) 
    550             END_2D 
    551          ENDIF    
     564               zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   & 
     565                  &                                      + e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) ) * ssumask(ji,jj) 
     566               zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   & 
     567                  &                                      + e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) ) * ssvmask(ji,jj) 
     568            END_2D 
     569#endif 
     570         ENDIF 
    552571         !          
    553572         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 
    … …  
    624643               !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 
    625644               !                    ! backward interpolated depth used in spg terms at jn+1/2 
     645#if defined key_qcoTest_FluxForm 
     646            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     647               zhu_bck = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     648               zhv_bck = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     649#else 
    626650               zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
    627651                    &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
    628652               zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
    629653                    &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     654#endif 
    630655               !                    ! inverse depth at jn+1 
    631656               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
    … …  
    646671         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 
    647672            DO_2D( 0, 0, 0, 0 ) 
    648                   ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 
    649                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 
     673               ua_e(ji,jj) =  ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) ) 
     674               va_e(ji,jj) =  va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) ) 
    650675            END_2D 
    651676         ENDIF 
    652677        
    653678         IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 
    654             hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
    655             hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
    656             hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
    657             hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
     679            hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
     680            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
     681            hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
     682            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
    658683         ENDIF 
    659684         ! 
    … …  
    743768      ELSE 
    744769         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 
     770#if defined key_qcoTest_FluxForm 
     771         !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
    745772         DO_2D( 1, 0, 1, 0 ) 
    746             zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) & 
    747                &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      & 
    748                &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 
    749             zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) & 
    750                &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      & 
    751                &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 
    752          END_2D 
     773            zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj) 
     774            zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj) 
     775         END_2D 
     776#else 
     777         DO_2D( 1, 0, 1, 0 ) 
     778            zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   & 
     779               &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 
     780            zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   & 
     781               &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 
     782         END_2D 
     783#endif    
    753784         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    754785         ! 
    755786         DO jk=1,jpkm1 
    756             puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 
    757             pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 
     787            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   & 
     788               &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 
     789            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   & 
     790               &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 
    758791         END DO 
    759792         ! Save barotropic velocities not transport: 
    … …  
    899932      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    900933         !                                   ! --------------- 
    901          IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file 
     934         IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN    !* Read the restart file 
    902935            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp )    
    903936            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp )  
    … …  
    10601093      !! although they should be updated in the variable volume case. Not a big approximation. 
    10611094      !! To remove this approximation, copy lines below inside barotropic loop 
    1062       !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 
     1095      !! and update depths at T- points (ht) at each barotropic time step 
    10631096      !! 
    10641097      !! Compute zwz = f / ( height of the water colomn ) 
    … …  
    10671100      INTEGER  ::   ji ,jj, jk              ! dummy loop indices 
    10681101      REAL(wp) ::   z1_ht 
    1069       REAL(wp), DIMENSION(jpi,jpj) :: zhf 
    10701102      !!---------------------------------------------------------------------- 
    10711103      ! 
    10721104      SELECT CASE( nvor_scheme ) 
    1073       CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme 
    1074          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
     1105      CASE( np_EEN, np_ENE, np_ENS , np_MIX )   !=  schemes using the same e3f definition 
     1106         SELECT CASE( nn_e3f_typ )                  !* ff_f/e3 at F-point 
    10751107         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    1076             DO_2D( 1, 0, 1, 0 ) 
    1077                zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
    1078                     &           ht(ji  ,jj  ) + ht(ji+1,jj  )  ) * 0.25_wp   
     1108            DO_2D( 0, 0, 0, 0 ) 
     1109               zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   & 
     1110                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp   
    10791111               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    10801112            END_2D 
    10811113         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
    1082             DO_2D( 1, 0, 1, 0 ) 
    1083                zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      & 
    1084                     &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   & 
    1085                     &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
    1086                     &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     1114            DO_2D( 0, 0, 0, 0 ) 
     1115               zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      & 
     1116                    &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   & 
     1117                    &    / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1)      & 
     1118                    &          + ssmask(ji,jj  ) + ssmask(ji+1,jj  ) , 1._wp  )   ) 
    10871119               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    10881120            END_2D 
    10891121         END SELECT 
    10901122         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
    1091          ! 
    1092          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1123      END SELECT 
     1124      ! 
     1125      SELECT CASE( nvor_scheme ) 
     1126      CASE( np_EEN ) 
     1127         ! 
     1128         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
    10931129         DO_2D( 0, 1, 0, 1 ) 
    10941130            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    … …  
    10981134         END_2D 
    10991135         ! 
    1100       CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme 
    1101          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1136      CASE( np_EET )                            != EEN scheme using e3t energy conserving scheme 
     1137         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;  ftsw(1,:) = 0._wp 
    11021138         DO_2D( 0, 1, 0, 1 ) 
    11031139            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 
    … …  
    11081144         END_2D 
    11091145         ! 
    1110       CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
    1111          ! 
    1112          zwz(:,:) = 0._wp 
    1113          zhf(:,:) = 0._wp 
    1114           
    1115          !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed  
    1116 !!gm    A priori a better value should be something like : 
    1117 !!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
    1118 !!gm                     divided by the sum of the corresponding mask  
    1119 !!gm  
    1120 !!             
    1121          IF( .NOT.ln_sco ) THEN 
    1122    
    1123    !!gm  agree the JC comment  : this should be done in a much clear way 
    1124    
    1125    ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    1126    !     Set it to zero for the time being  
    1127    !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
    1128    !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
    1129    !              ENDIF 
    1130    !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    1131             ! 
    1132          ELSE 
    1133             ! 
    1134             !zhf(:,:) = hbatf(:,:) 
    1135             DO_2D( 1, 0, 1, 0 ) 
    1136                zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
    1137                     &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
    1138                     &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
    1139                     &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
    1140             END_2D 
    1141          ENDIF 
    1142          ! 
    1143          DO jj = 1, jpjm1 
    1144             zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    1145          END DO 
    1146          ! 
    1147          DO jk = 1, jpkm1 
    1148             DO jj = 1, jpjm1 
    1149                zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
    1150             END DO 
    1151          END DO 
    1152          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 
    1153          ! JC: TBC. hf should be greater than 0  
    1154          DO_2D( 1, 1, 1, 1 ) 
    1155             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) 
    1156          END_2D 
    1157          zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    11581146      END SELECT 
    11591147       
    11601148   END SUBROUTINE dyn_cor_2d_init 
    1161  
    11621149 
    11631150 
    … …  
    13531340   END SUBROUTINE wad_spg 
    13541341      
    1355  
    13561342 
    13571343   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/dynvor.F90

    r14049 r14054  
    2121   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 
    2222   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation 
     23   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity 
    2324   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T) 
    2425   !!---------------------------------------------------------------------- 
    … …  
    2627   !!---------------------------------------------------------------------- 
    2728   !!   dyn_vor       : Update the momentum trend with the vorticity trend 
     29   !!       vor_enT   : energy conserving scheme at T-pt  (ln_dynvor_enT=T) 
     30   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T) 
    2831   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T) 
    29    !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T) 
    3032   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T) 
     33   !!       vor_eeT   : energy conserving at T-pt         (ln_dynvor_eeT=T) 
    3134   !!   dyn_vor_init  : set and control of the different vorticity option 
    3235   !!---------------------------------------------------------------------- 
    … …  
    5861   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
    5962   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    60    INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    6163   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX) 
    6264   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     65   INTEGER, PUBLIC ::   nn_e3f_typ      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    6366 
    6467   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme 
    … …  
    8184   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation 
    8285   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
    83    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation 
    84    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       -  
     86   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
     87   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
     88   ! 
     89   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   e3f_0vor   ! e3f used in EEN, ENE and ENS cases (key_qco only) 
    8590    
    8691   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
    … …  
    235240      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    236241      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
    237       REAL(wp), DIMENSION(jpi,jpj)       ::   zwx, zwy, zwt   ! 2D workspace 
    238       REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz             ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 
     242      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace 
     243      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwz      ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 
    239244      !!---------------------------------------------------------------------- 
    240245      ! 
    … …  
    246251      ! 
    247252      ! 
    248       SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
    249       CASE ( np_RVO )                           !* relative vorticity 
    250          DO jk = 1, jpkm1                                 ! Horizontal slab 
     253      SELECT CASE( kvor )                 !== relative vorticity considered  ==! 
     254      ! 
     255      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used 
     256         ALLOCATE( zwz(jpi,jpj,jpk) ) 
     257         DO jk = 1, jpkm1                                ! Horizontal slab 
    251258            DO_2D( 1, 0, 1, 0 ) 
    252259               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    253260                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    254261            END_2D 
    255             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     262            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity  
    256263               DO_2D( 1, 0, 1, 0 ) 
    257264                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    … …  
    259266            ENDIF 
    260267         END DO 
    261  
    262268         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    263  
    264       CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    265          DO jk = 1, jpkm1                                 ! Horizontal slab 
    266             DO_2D( 1, 0, 1, 0 )                          ! relative vorticity 
    267                zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
    268                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    269             END_2D 
    270             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
    271                DO_2D( 1, 0, 1, 0 ) 
    272                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    273                END_2D 
    274             ENDIF 
    275          END DO 
    276  
    277          CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    278  
     269         ! 
    279270      END SELECT 
    280271 
    281272      !                                                ! =============== 
    282273      DO jk = 1, jpkm1                                 ! Horizontal slab 
    283       !                                                ! =============== 
    284  
     274         !                                             ! =============== 
     275         ! 
    285276         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
     277         ! 
    286278         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    287279            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 
    288280         CASE ( np_RVO )                           !* relative vorticity 
    289281            DO_2D( 0, 1, 0, 1 ) 
    290                zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    291                   &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 
    292                   &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     282               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       & 
     283                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )  & 
     284                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    293285            END_2D 
    294286         CASE ( np_MET )                           !* metric term 
    295287            DO_2D( 0, 1, 0, 1 ) 
    296                zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    297                   &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) & 
    298                   &             * e3t(ji,jj,jk,Kmm) 
     288               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       & 
     289                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   & 
     290                  &       * e3t(ji,jj,jk,Kmm) 
    299291            END_2D 
    300292         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    301293            DO_2D( 0, 1, 0, 1 ) 
    302                zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    303                   &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) & 
    304                   &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     294               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        & 
     295                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   & 
     296                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    305297            END_2D 
    306298         CASE ( np_CME )                           !* Coriolis + metric 
    307299            DO_2D( 0, 1, 0, 1 ) 
    308                zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    309                     &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    310                     &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) & 
    311                     &          * e3t(ji,jj,jk,Kmm) 
     300               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               & 
     301                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      & 
     302                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   & 
     303                    &     * e3t(ji,jj,jk,Kmm) 
    312304            END_2D 
    313305         CASE DEFAULT                                             ! error 
    314             CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     306            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor') 
    315307         END SELECT 
    316308         ! 
    … …  
    328320      END DO                                           !   End of slab 
    329321      !                                                ! =============== 
     322      ! 
     323      SELECT CASE( kvor )        ! deallocate zwz if necessary 
     324      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz ) 
     325      END SELECT 
     326      ! 
    330327   END SUBROUTINE vor_enT 
    331328 
    … …  
    358355      ! 
    359356      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    360       REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
     357      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars 
    361358      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace 
    362359      !!---------------------------------------------------------------------- 
    … …  
    380377                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    381378            END_2D 
     379            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     380               DO_2D( 1, 0, 1, 0 ) 
     381                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     382               END_2D 
     383            ENDIF 
    382384         CASE ( np_MET )                           !* metric term 
    383385            DO_2D( 1, 0, 1, 0 ) 
    … …  
    390392                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    391393            END_2D 
     394            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term) 
     395               DO_2D( 1, 0, 1, 0 ) 
     396                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
     397               END_2D 
     398            ENDIF 
    392399         CASE ( np_CME )                           !* Coriolis + metric 
    393400            DO_2D( 1, 0, 1, 0 ) 
    … …  
    399406         END SELECT 
    400407         ! 
    401          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    402             DO_2D( 1, 0, 1, 0 ) 
    403                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    404             END_2D 
    405          ENDIF 
    406  
    407          IF( ln_sco ) THEN 
    408             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    409             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    410             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    411          ELSE 
    412             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    413             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    414          ENDIF 
     408#if defined key_qco 
     409         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco) 
     410            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk) 
     411         END_2D 
     412#else 
     413         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==! 
     414         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     415            DO_2D( 1, 0, 1, 0 ) 
     416               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     417                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     418                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     419                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     420               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f 
     421               ELSE                       ;   zwz(ji,jj) = 0._wp 
     422               ENDIF 
     423            END_2D 
     424         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     425            DO_2D( 1, 0, 1, 0 ) 
     426               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     427                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     428                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     429                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   ) 
     430               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   & 
     431                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   ) 
     432               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f 
     433               ELSE                       ;   zwz(ji,jj) = 0._wp 
     434               ENDIF 
     435            END_2D 
     436         END SELECT 
     437#endif 
     438         !                                   !==  horizontal fluxes  ==! 
     439         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     440         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     441         ! 
    415442         !                                   !==  compute and add the vorticity term trend  =! 
    416443         DO_2D( 0, 0, 0, 0 ) 
    … …  
    455482      ! 
    456483      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    457       REAL(wp) ::   zuav, zvau   ! local scalars 
     484      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars 
    458485      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace 
    459486      !!---------------------------------------------------------------------- 
    … …  
    476503                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    477504            END_2D 
     505            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     506               DO_2D( 1, 0, 1, 0 ) 
     507                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     508               END_2D 
     509            ENDIF 
    478510         CASE ( np_MET )                           !* metric term 
    479511            DO_2D( 1, 0, 1, 0 ) 
    … …  
    486518                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    487519            END_2D 
     520            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term) 
     521               DO_2D( 1, 0, 1, 0 ) 
     522                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
     523               END_2D 
     524            ENDIF 
    488525         CASE ( np_CME )                           !* Coriolis + metric 
    489526            DO_2D( 1, 0, 1, 0 ) 
    … …  
    495532         END SELECT 
    496533         ! 
    497          IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==! 
    498             DO_2D( 1, 0, 1, 0 ) 
    499                zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    500             END_2D 
    501          ENDIF 
    502          ! 
    503          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    504             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    505             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    506             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    507          ELSE 
    508             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    509             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    510          ENDIF 
     534         ! 
     535#if defined key_qco 
     536         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco) 
     537            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk) 
     538         END_2D 
     539#else 
     540         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==! 
     541         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     542            DO_2D( 1, 0, 1, 0 ) 
     543               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     544                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     545                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     546                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     547               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f 
     548               ELSE                       ;   zwz(ji,jj) = 0._wp 
     549               ENDIF 
     550            END_2D 
     551         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     552            DO_2D( 1, 0, 1, 0 ) 
     553               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     554                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     555                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     556                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   ) 
     557               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   & 
     558                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   ) 
     559               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f 
     560               ELSE                       ;   zwz(ji,jj) = 0._wp 
     561               ENDIF 
     562            END_2D 
     563         END SELECT 
     564#endif 
     565         !                                   !==  horizontal fluxes  ==! 
     566         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     567         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     568         ! 
    511569         !                                   !==  compute and add the vorticity term trend  =! 
    512570         DO_2D( 0, 0, 0, 0 ) 
    … …  
    566624         !                                             ! =============== 
    567625         ! 
    568          SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point 
     626#if defined key_qco 
     627         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco) 
     628            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk) 
     629         END_2D 
     630#else 
     631         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point 
    569632         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    570633            DO_2D( 1, 0, 1, 0 ) 
    … …  
    590653            END_2D 
    591654         END SELECT 
     655#endif 
    592656         ! 
    593657         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     658         ! 
    594659         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    595660            DO_2D( 1, 0, 1, 0 ) 
    … …  
    601666                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 
    602667            END_2D 
     668            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     669               DO_2D( 1, 0, 1, 0 ) 
     670                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     671               END_2D 
     672            ENDIF 
    603673         CASE ( np_MET )                           !* metric term 
    604674            DO_2D( 1, 0, 1, 0 ) 
    … …  
    612682                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    613683            END_2D 
     684            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     685               DO_2D( 1, 0, 1, 0 ) 
     686                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)  
     687               END_2D 
     688            ENDIF 
    614689         CASE ( np_CME )                           !* Coriolis + metric 
    615690            DO_2D( 1, 0, 1, 0 ) 
    … …  
    620695            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    621696         END SELECT 
    622          ! 
    623          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    624             DO_2D( 1, 0, 1, 0 ) 
    625                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    626             END_2D 
    627          ENDIF 
     697         !                                             ! =============== 
    628698      END DO                                           !   End of slab 
    629          ! 
     699      !                                                ! =============== 
     700      ! 
    630701      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    631  
     702      ! 
     703      !                                                ! =============== 
    632704      DO jk = 1, jpkm1                                 ! Horizontal slab 
     705         !                                             ! =============== 
    633706         ! 
    634707         !                                   !==  horizontal fluxes  ==! 
    635708         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    636709         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    637  
     710         ! 
    638711         !                                   !==  compute and add the vorticity term trend  =! 
    639          jj = 2 
    640          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    641          DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    642                ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
    643                ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
    644                ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
    645                ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
    646          END DO 
    647          DO jj = 3, jpj 
    648             DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3 
    649                ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
    650                ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
    651                ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
    652                ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
    653             END DO 
    654          END DO 
     712         DO_2D( 0, 1, 0, 1 ) 
     713            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) 
     714            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) 
     715            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) 
     716            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) 
     717         END_2D 
     718         ! 
    655719         DO_2D( 0, 0, 0, 0 ) 
    656720            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    … …  
    667731 
    668732 
    669  
    670733   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs ) 
    671734      !!---------------------------------------------------------------------- 
    … …  
    685748      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    686749      !!---------------------------------------------------------------------- 
    687       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     750      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index 
    688751      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    689       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    690       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities 
    691       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
     752      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric 
     753      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities 
     754      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend 
    692755      ! 
    693756      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    … …  
    702765      IF( kt == nit000 ) THEN 
    703766         IF(lwp) WRITE(numout,*) 
    704          IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
     767         IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme' 
    705768         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    706769      ENDIF 
    … …  
    722785                  &          * r1_e1e2f(ji,jj) 
    723786            END_2D 
     787            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     788               DO_2D( 1, 0, 1, 0 ) 
     789                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     790               END_2D 
     791            ENDIF 
    724792         CASE ( np_MET )                           !* metric term 
    725793            DO_2D( 1, 0, 1, 0 ) 
    … …  
    733801                  &                         * r1_e1e2f(ji,jj)    ) 
    734802            END_2D 
     803            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
     804               DO_2D( 1, 0, 1, 0 ) 
     805                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)  
     806               END_2D 
     807            ENDIF 
    735808         CASE ( np_CME )                           !* Coriolis + metric 
    736809            DO_2D( 1, 0, 1, 0 ) 
    … …  
    742815         END SELECT 
    743816         ! 
    744          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    745             DO_2D( 1, 0, 1, 0 ) 
    746                zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    747             END_2D 
    748          ENDIF 
    749       END DO 
     817         !                                             ! =============== 
     818      END DO                                           !   End of slab 
     819      !                                                ! =============== 
    750820      ! 
    751821      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 
    752822      ! 
     823      !                                                ! =============== 
    753824      DO jk = 1, jpkm1                                 ! Horizontal slab 
    754  
    755       !                                   !==  horizontal fluxes  ==! 
     825         !                                             ! =============== 
     826         ! 
     827         !                                   !==  horizontal fluxes  ==! 
    756828         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    757829         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    758  
     830         ! 
    759831         !                                   !==  compute and add the vorticity term trend  =! 
    760          jj = 2 
    761          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
    762          DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
    763                z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    764                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    765                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    766                ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
    767                ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
    768          END DO 
    769          DO jj = 3, jpj 
    770             DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3 
    771                z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
    772                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
    773                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
    774                ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
    775                ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
    776             END DO 
    777          END DO 
     832         DO_2D( 0, 1, 0, 1 ) 
     833            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
     834            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t 
     835            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t 
     836            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t 
     837            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t 
     838         END_2D 
     839         ! 
    778840         DO_2D( 0, 0, 0, 0 ) 
    779841            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    … …  
    799861      INTEGER ::   ji, jj, jk    ! dummy loop indices 
    800862      INTEGER ::   ioptio, ios   ! local integer 
     863      REAL(wp) ::   zmsk    ! local scalars 
    801864      !! 
    802865      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
    803          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     866         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk 
    804867      !!---------------------------------------------------------------------- 
    805868      ! 
    … …  
    823886         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT 
    824887         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
    825          WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     888         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ 
    826889         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
    827890         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    828891      ENDIF 
    829  
    830       IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available') 
    831892 
    832893!!gm  this should be removed when choosing a unique strategy for fmask at the coast 
    … …  
    891952         ! 
    892953      END SELECT 
    893        
     954#if defined key_qco 
     955      SELECT CASE( nvor_scheme )    ! qco case: pre-computed a specific e3f_0 for some vorticity schemes 
     956      CASE( np_ENS , np_ENE , np_EEN , np_MIX ) 
     957         ! 
     958         ALLOCATE( e3f_0vor(jpi,jpj,jpk) ) 
     959         ! 
     960         SELECT CASE( nn_e3f_typ ) 
     961         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4) 
     962            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     963               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   & 
     964                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     965                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   & 
     966                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp 
     967            END_3D 
     968         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     969            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     970               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   & 
     971                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  ) 
     972               ! 
     973               IF( zmsk /= 0._wp ) THEN  
     974                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   & 
     975                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     976                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   & 
     977                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk 
     978               ENDIF 
     979            END_3D 
     980         END SELECT 
     981         ! 
     982         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp ) 
     983         !                                 ! insure e3f_0vor /= 0 
     984         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:) 
     985         ! 
     986      END SELECT 
     987      ! 
     988#endif 
    894989      IF(lwp) THEN                   ! Print the choice 
    895990         WRITE(numout,*) 
    … …  
    898993         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
    899994         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
     995                              IF( ln_dynadv_vec )   CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form') 
    900996         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
    901997         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/DYN/sshwzv.F90

    r13497 r14054  
    66   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code 
    77   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA  
    8    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    9    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    10    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
    11    !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection 
    12    !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 
     8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea)  Assimilation interface 
     9   !!             -   !  2010-09  (D.Storkey and E.O'Dea)  bug fixes for BDY module 
     10   !!            3.3  !  2011-10  (M. Leclair)  split former ssh_wzv routine and remove all vvl related work 
     11   !!            4.0  !  2018-12  (A. Coward)  add mixed implicit/explicit advection 
     12   !!            4.1  !  2019-08  (A. Coward, D. Storkey)  Rename ssh_nxt -> ssh_atf. Now only does time filtering. 
     13   !!             -   !  2020-08  (S. Techene, G. Madec)  add here ssh initiatlisation 
    1314   !!---------------------------------------------------------------------- 
    1415 
    … …  
    1718   !!   ssh_atf       : time filter the ssh arrays 
    1819   !!   wzv           : compute now vertical velocity 
     20   !!   ssh_init_rst  : ssh set from restart or domcfg.nc file or usr_def_istat_ssh 
    1921   !!---------------------------------------------------------------------- 
    2022   USE oce            ! ocean dynamics and tracers variables 
    … …  
    4042   USE timing         ! Timing 
    4143   USE wet_dry        ! Wetting/Drying flux limiting 
    42  
     44   USE usrdef_istate, ONLY : usr_def_istate_ssh   ! user defined ssh initial state  
     45    
    4346   IMPLICIT NONE 
    4447   PRIVATE 
    4548 
    46    PUBLIC   ssh_nxt    ! called by step.F90 
    47    PUBLIC   wzv        ! called by step.F90 
    48    PUBLIC   wAimp      ! called by step.F90 
    49    PUBLIC   ssh_atf    ! called by step.F90 
     49   PUBLIC   ssh_nxt        ! called by step.F90 
     50   PUBLIC   wzv            ! called by step.F90 
     51   PUBLIC   wAimp          ! called by step.F90 
     52   PUBLIC   ssh_atf        ! called by step.F90 
     53   PUBLIC   ssh_init_rst   ! called by domain.F90 
    5054 
    5155   !! * Substitutions 
    5256#  include "do_loop_substitute.h90" 
    5357#  include "domzgr_substitute.h90" 
    54  
    5558   !!---------------------------------------------------------------------- 
    5659   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    299302         !                                                  ! filtered "now" field 
    300303         pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 
     304         ! 
    301305         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed 
    302306            zcoef = rn_atfp * rn_Dt * r1_rho0 
    … …  
    307311 
    308312            ! ice sheet coupling 
    309             IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 
     313            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   & 
     314               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 
    310315 
    311316         ENDIF 
    312317      ENDIF 
    313318      ! 
    314       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask ) 
     319      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask ) 
    315320      ! 
    316321      IF( ln_timing )   CALL timing_stop('ssh_atf') 
    … …  
    431436      ! 
    432437   END SUBROUTINE wAimp 
     438 
     439 
     440   SUBROUTINE ssh_init_rst( Kbb, Kmm, Kaa ) 
     441      !!--------------------------------------------------------------------- 
     442      !!                   ***  ROUTINE ssh_init_rst  *** 
     443      !! 
     444      !! ** Purpose :   ssh initialization of the sea surface height (ssh) 
     445      !! 
     446      !! ** Method  :   set ssh from restart or read configuration, or user_def 
     447      !!              * ln_rstart = T 
     448      !!                   USE of IOM library to read ssh in the restart file 
     449      !!                   Leap-Frog: Kbb and Kmm are read except for l_1st_euler=T 
     450      !! 
     451      !!              * otherwise  
     452      !!                   call user defined ssh or 
     453      !!                   set to -ssh_ref in wet and drying case with domcfg.nc 
     454      !! 
     455      !!              NB: ssh_b/n are written by restart.F90 
     456      !!---------------------------------------------------------------------- 
     457      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time level indices 
     458      ! 
     459      INTEGER ::   ji, jj, jk 
     460      !!---------------------------------------------------------------------- 
     461      ! 
     462      IF(lwp) THEN 
     463         WRITE(numout,*) 
     464         WRITE(numout,*) 'ssh_init_rst : ssh initialization' 
     465         WRITE(numout,*) '~~~~~~~~~~~~ ' 
     466      ENDIF 
     467      ! 
     468      !                            !=============================! 
     469      IF( ln_rstart ) THEN         !==  Read the restart file  ==! 
     470         !                         !=============================! 
     471         ! 
     472         !                                     !*  Read ssh at Kmm 
     473         IF(lwp) WRITE(numout,*) 
     474         IF(lwp) WRITE(numout,*)    '      Kmm sea surface height read in the restart file' 
     475         CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm) ) 
     476         ! 
     477         IF( l_1st_euler ) THEN                !* Euler at first time-step 
     478            IF(lwp) WRITE(numout,*) 
     479            IF(lwp) WRITE(numout,*) '      Euler first time step : ssh(Kbb) = ssh(Kmm)' 
     480            ssh(:,:,Kbb) = ssh(:,:,Kmm) 
     481            ! 
     482         ELSE                                  !* read ssh at Kbb 
     483            IF(lwp) WRITE(numout,*) 
     484            IF(lwp) WRITE(numout,*) '      Kbb sea surface height read in the restart file' 
     485            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) ) 
     486         ENDIF 
     487         !                         !============================! 
     488      ELSE                         !==  Initialize at "rest"  ==! 
     489         !                         !============================! 
     490         ! 
     491         IF(lwp) WRITE(numout,*) 
     492         IF(lwp) WRITE(numout,*)    '      initialization at rest' 
     493         ! 
     494         IF( ll_wd ) THEN                      !* wet and dry  
     495            ! 
     496            IF( ln_read_cfg  ) THEN                 ! read configuration : ssh_ref is read in domain_cfg file 
     497!!st  why ssh is not masked : i.e. ssh(:,:,Kmm) = -ssh_ref*ssmask(:,:), 
     498!!st  since at the 1st time step lbclnk will be applied on ssh at Kaa but not initially at Kbb and Kmm 
     499               ssh(:,:,Kbb) = -ssh_ref 
     500               ! 
     501               DO_2D( 1, 1, 1, 1 ) 
     502                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN   ! if total depth is less than min depth 
     503                     ssh(ji,jj,Kbb) = rn_wdmin1 - ht_0(ji,jj) 
     504                  ENDIF 
     505               END_2D 
     506            ELSE                                    ! user define configuration case   
     507               CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) ) 
     508            ENDIF 
     509            ! 
     510         ELSE                                  !* user defined configuration 
     511            CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) ) 
     512            ! 
     513         ENDIF 
     514         ! 
     515         ssh(:,:,Kmm) = ssh(:,:,Kbb)           !* set now values from to before ones 
     516         ssh(:,:,Kaa) = 0._wp  
     517      ENDIF 
     518      ! 
     519   END SUBROUTINE ssh_init_rst 
     520       
    433521   !!====================================================================== 
    434522END MODULE sshwzv 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icb_oce.F90

    r13286 r14054  
    5757   TYPE, PUBLIC ::   point              !: properties of an individual iceberg (position, mass, size, etc...) 
    5858      INTEGER  ::   year 
    59       REAL(wp) ::   xi , yj                                                   ! iceberg coordinates in the (i,j) referential (global) 
     59      REAL(wp) ::   xi , yj , zk                                              ! iceberg coordinates in the (i,j) referential (global) and deepest level affected 
    6060      REAL(wp) ::   e1 , e2                                                   ! horizontal scale factors at the iceberg position 
    6161      REAL(wp) ::   lon, lat, day                                             ! geographic position 
    6262      REAL(wp) ::   mass, thickness, width, length, uvel, vvel                ! iceberg physical properties 
    63       REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi, sss    ! properties of iceberg environment  
     63      REAL(wp) ::   ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, sss, cn, hi  ! properties of iceberg environment  
    6464      REAL(wp) ::   mass_of_bits, heat_density 
     65      INTEGER  ::   kb                                                   ! icb bottom level 
    6566   END TYPE point 
    6667 
    … …  
    8586   ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position 
    8687   ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 
    87    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e 
    88    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e, ss_e 
     88   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssu_e, ssv_e 
     89   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   sst_e, sss_e, fr_e 
    8990   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e 
    9091   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e 
    9192   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tmask_e, umask_e, vmask_e 
     93   REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   rlon_e, rlat_e, ff_e 
     94   REAl(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   uoce_e, voce_e, toce_e, e3t_e 
     95   ! 
    9296#if defined key_si3 || defined key_cice 
    9397   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   hi_e, ui_e, vi_e 
    … …  
    117121   INTEGER , PUBLIC ::   nn_verbose_write                !: timesteps between verbose messages 
    118122   REAL(wp), PUBLIC ::   rn_rho_bergs                    !: Density of icebergs 
     123   REAL(wp), PUBLIC ::   rho_berg_1_oce                  !: convertion factor (thickness to draft) (rn_rho_bergs/pp_rho_seawater) 
    119124   REAL(wp), PUBLIC ::   rn_LoW_ratio                    !: Initial ratio L/W for newly calved icebergs 
    120125   REAL(wp), PUBLIC ::   rn_bits_erosion_fraction        !: Fraction of erosion melt flux to divert to bergy bits 
    … …  
    124129   LOGICAL , PUBLIC ::   ln_time_average_weight          !: Time average the weight on the ocean    !!gm I don't understand that ! 
    125130   REAL(wp), PUBLIC ::   rn_speed_limit                  !: CFL speed limit for a berg 
     131   LOGICAL , PUBLIC ::   ln_M2016, ln_icb_grd            !: use Nacho's Merino 2016 work 
    126132   ! 
    127133   ! restart 
    … …  
    135141   REAL(wp), DIMENSION(nclasses), PUBLIC ::   rn_initial_thickness !  Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run 
    136142   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   src_calving, src_calving_hflx    !: accumulate input ice 
     143   INTEGER , PUBLIC             , SAVE                     ::   micbkb                           !: deepest level affected by icebergs 
    137144   INTEGER , PUBLIC             , SAVE                     ::   numicb                           !: iceberg IO 
    138145   INTEGER , PUBLIC             , SAVE, DIMENSION(nkounts) ::   num_bergs                        !: iceberg counter 
    … …  
    171178      ! 
    172179      ! expanded arrays for bilinear interpolation 
    173       ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,    & 
    174          &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,    & 
     180      ALLOCATE( ssu_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   & 
     181         &      ssv_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   & 
    175182#if defined key_si3 || defined key_cice 
    176183         &      ui_e(0:jpi+1,0:jpj+1) ,                            & 
    … …  
    178185         &      hi_e(0:jpi+1,0:jpj+1) ,                            & 
    179186#endif 
    180          &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   & 
    181          &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   & 
    182          &      ss_e(0:jpi+1,0:jpj+1) ,                            &  
     187         &      fr_e(0:jpi+1,0:jpj+1) ,                            & 
     188         &      sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,  & 
     189         &      sss_e(0:jpi+1,0:jpj+1) ,                           &  
    183190         &      first_width(nclasses) , first_length(nclasses) ,   & 
    184191         &      src_calving (jpi,jpj) ,                            & 
    … …  
    186193      icb_alloc = icb_alloc + ill 
    187194 
     195      IF ( ln_M2016 ) THEN 
     196         ALLOCATE( uoce_e(0:jpi+1,0:jpj+1,jpk), voce_e(0:jpi+1,0:jpj+1,jpk), & 
     197            &      toce_e(0:jpi+1,0:jpj+1,jpk), e3t_e(0:jpi+1,0:jpj+1,jpk) , STAT=ill ) 
     198         icb_alloc = icb_alloc + ill 
     199      END IF 
     200      ! 
    188201      ALLOCATE( tmask_e(0:jpi+1,0:jpj+1), umask_e(0:jpi+1,0:jpj+1), vmask_e(0:jpi+1,0:jpj+1), & 
    189          &      STAT=ill) 
     202         &      rlon_e(0:jpi+1,0:jpj+1) , rlat_e(0:jpi+1,0:jpj+1) , ff_e(0:jpi+1,0:jpj+1)   , STAT=ill) 
    190203      icb_alloc = icb_alloc + ill 
    191204 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icbclv.F90

    r13295 r14054  
    2121   USE lbclnk         ! NEMO boundary exchanges for gridded data 
    2222 
    23    USE icb_oce        ! iceberg variables 
    2423   USE icbdia         ! iceberg diagnostics 
    2524   USE icbutl         ! iceberg utility routines 
    … …  
    142141                  newpt%mass           = rn_initial_mass     (jn) 
    143142                  newpt%thickness      = rn_initial_thickness(jn) 
     143                  newpt%kb             = 1         ! compute correctly in icbthm if needed         
    144144                  newpt%width          = first_width         (jn) 
    145145                  newpt%length         = first_length        (jn) 
    … …  
    172172      END DO 
    173173      ! 
    174       DO jn = 1, nclasses 
    175          CALL lbc_lnk( 'icbclv', berg_grid%stored_ice(:,:,jn), 'T', 1._wp ) 
    176       END DO 
    177       CALL lbc_lnk( 'icbclv', berg_grid%stored_heat, 'T', 1._wp ) 
    178       ! 
    179174      IF( nn_verbose_level > 0 .AND. icntmax > 1 )   WRITE(numicb,*) 'icb_clv: icnt=', icnt,' on', narea 
    180175      ! 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icbdyn.F90

    r13281 r14054  
    1414   USE dom_oce        ! NEMO ocean domain 
    1515   USE phycst         ! NEMO physical constants 
     16   USE in_out_manager                      ! IO parameters 
    1617   ! 
    1718   USE icb_oce        ! define iceberg arrays 
    … …  
    9798         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1 
    9899         ! 
    99          CALL icb_ground( zxi2, zxi1, zu1,   & 
    100             &             zyj2, zyj1, zv1, ll_bounced ) 
     100         CALL icb_ground( berg, zxi2, zxi1, zu1,   & 
     101            &                   zyj2, zyj1, zv1, ll_bounced ) 
    101102 
    102103         !                                         !**   A2 = A(X2,V2) 
    … …  
    113114         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
    114115         ! 
    115          CALL icb_ground( zxi3, zxi1, zu3,   & 
    116             &                zyj3, zyj1, zv3, ll_bounced ) 
     116         CALL icb_ground( berg, zxi3, zxi1, zu3,   & 
     117            &                   zyj3, zyj1, zv3, ll_bounced ) 
    117118 
    118119         !                                         !**   A3 = A(X3,V3) 
    … …  
    129130         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
    130131 
    131          CALL icb_ground( zxi4, zxi1, zu4,   & 
    132             &             zyj4, zyj1, zv4, ll_bounced ) 
     132         CALL icb_ground( berg, zxi4, zxi1, zu4,   & 
     133            &                   zyj4, zyj1, zv4, ll_bounced ) 
    133134 
    134135         !                                         !**   A4 = A(X4,V4) 
    … …  
    148149         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 ) 
    149150 
    150          CALL icb_ground( zxi_n, zxi1, zuvel_n,   & 
    151             &             zyj_n, zyj1, zvvel_n, ll_bounced ) 
     151         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   & 
     152            &                   zyj_n, zyj1, zvvel_n, ll_bounced ) 
    152153 
    153154         pt%uvel = zuvel_n                        !** save in berg structure 
    … …  
    156157         pt%yj   = zyj_n 
    157158 
    158          ! update actual position 
    159          pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 
    160          pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 
    161  
    162159         berg => berg%next                         ! switch to the next berg 
    163160         ! 
    … …  
    167164 
    168165 
    169    SUBROUTINE icb_ground( pi, pi0, pu,   & 
    170       &                   pj, pj0, pv, ld_bounced ) 
     166   SUBROUTINE icb_ground( berg, pi, pi0, pu,   & 
     167      &                         pj, pj0, pv, ld_bounced ) 
    171168      !!---------------------------------------------------------------------- 
    172169      !!                  ***  ROUTINE icb_ground  *** 
    … …  
    177174      !!                NB two possibilities available one of which is hard-coded here 
    178175      !!---------------------------------------------------------------------- 
     176      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg 
     177      ! 
    179178      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position 
    180179      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position 
    … …  
    184183      INTEGER  ::   ii, ii0 
    185184      INTEGER  ::   ij, ij0 
     185      INTEGER  ::   ikb 
    186186      INTEGER  ::   ibounce_method 
     187      ! 
     188      REAL(wp) :: zD  
     189      REAL(wp), DIMENSION(jpk) :: ze3t 
    187190      !!---------------------------------------------------------------------- 
    188191      ! 
    … …  
    200203      ij  = mj1( ij  ) 
    201204      ! 
    202       IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     205      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 
     206      IF ( ln_M2016 .AND. ln_icb_grd ) THEN 
     207         ! 
     208         ! draught (keel depth) 
     209         zD = rho_berg_1_oce * berg%current_point%thickness 
     210         ! 
     211         ! interpol needed data 
     212         CALL icb_utl_interp( pi, pj, pe3t=ze3t ) 
     213         !  
     214         !compute bottom level 
     215         CALL icb_utl_getkb( ikb, ze3t, zD ) 
     216         ! 
     217         ! berg reach a new t-cell, but an ocean one 
     218         ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0) 
     219         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN 
     220         ! 
     221      ELSE 
     222         IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     223      END IF 
    203224      ! 
    204225      ! From here, berg have reach land: treat grounding/bouncing 
    … …  
    257278      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    ! 
    258279      ! 
    259       INTEGER  ::   itloop 
    260       REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss 
    261       REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y 
     280      INTEGER  ::   itloop, ikb, jk 
     281      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 
     282      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y 
    262283      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF 
    263284      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 
    264       REAL(wp) ::   z_ocn, z_atm, z_ice 
     285      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep 
    265286      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 
    266287      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 
    267288      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 
     289      REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 
    268290      !!---------------------------------------------------------------------- 
    269291 
    270292      ! Interpolate gridded fields to berg 
    271293      nknberg = berg%number(1) 
    272       CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     & 
    273          &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss ) 
     294      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor 
     295         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities 
     296         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities 
     297         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient 
     298         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis 
    274299 
    275300      zM = berg%current_point%mass 
    276301      zT = berg%current_point%thickness               ! total thickness 
    277       zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth) 
     302      zD = rho_berg_1_oce * zT                        ! draught (keel depth) 
    278303      zF = zT - zD                                    ! freeboard 
    279304      zW = berg%current_point%width 
    … …  
    282307      zhi   = MIN( zhi   , zD    ) 
    283308      zD_hi = MAX( 0._wp, zD-zhi ) 
    284  
    285       ! Wave radiation 
    286       zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model 
     309  
     310     ! Wave radiation 
     311      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model 
    287312      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current; 
    288313      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here. 
    … …  
    309334      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp 
    310335 
     336      ! lateral velocities 
     337      ! default ssu and ssv 
     338      ! ln_M2016: mean velocity along the profile 
     339      IF ( ln_M2016 ) THEN 
     340         ! interpol needed data 
     341         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities 
     342         
     343         !compute bottom level 
     344         CALL icb_utl_getkb( ikb, ze3t, zD ) 
     345          
     346         ! compute mean velocity  
     347         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb) 
     348         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb) 
     349      ELSE 
     350         zuo = zssu 
     351         zvo = zssv 
     352      END IF 
     353 
    311354      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel 
    312355      ! 
    … …  
    321364         ! Explicit accelerations 
    322365         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 
    323          !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 
     366         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 
    324367         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 
    325          !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 
     368         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 
    326369         zaxe = -grav * zssh_x + zwave_rad * zuwave 
    327370         zaye = -grav * zssh_y + zwave_rad * zvwave 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icbini.F90

    r13295 r14054  
    7373      ! 
    7474      IF( .NOT. ln_icebergs )   RETURN 
    75  
     75      ! 
    7676      !                          ! allocate gridded fields 
    7777      IF( icb_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' ) 
    7878      ! 
    7979      !                          ! initialised variable with extra haloes to zero 
    80       uo_e(:,:) = 0._wp   ;   vo_e(:,:) = 0._wp   ; 
    81       ua_e(:,:) = 0._wp   ;   va_e(:,:) = 0._wp   ; 
    82       ff_e(:,:) = 0._wp   ;   tt_e(:,:) = 0._wp   ; 
    83       fr_e(:,:) = 0._wp   ;   ss_e(:,:) = 0._wp   ; 
     80      ssu_e(:,:) = 0._wp   ;   ssv_e(:,:) = 0._wp   ; 
     81      ua_e(:,:)  = 0._wp   ;   va_e(:,:)  = 0._wp   ; 
     82      ff_e(:,:)  = 0._wp   ;   sst_e(:,:) = 0._wp   ; 
     83      fr_e(:,:)  = 0._wp   ;   sss_e(:,:) = 0._wp   ; 
     84      ! 
     85      IF ( ln_M2016 ) THEN 
     86         toce_e(:,:,:) = 0._wp 
     87         uoce_e(:,:,:) = 0._wp 
     88         voce_e(:,:,:) = 0._wp 
     89         e3t_e(:,:,:)  = 0._wp 
     90      END IF 
     91      ! 
    8492#if defined key_si3 
    8593      hi_e(:,:) = 0._wp   ; 
    … …  
    100108      first_width (:) = SQRT(  rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) )  ) 
    101109      first_length(:) = rn_LoW_ratio * first_width(:) 
     110      rho_berg_1_oce  = rn_rho_bergs / pp_rho_seawater  ! scale factor used for convertion thickness to draft 
     111      ! 
     112      ! deepest level affected by icebergs 
     113      ! can be tuned but the safest is this  
     114      ! (with z* and z~ the depth of each level change overtime, so the more robust micbkb is jpk) 
     115      micbkb = jpk 
    102116 
    103117      berg_grid%calving      (:,:)   = 0._wp 
    … …  
    240254      vmask_e(:,:) = 0._wp   ;   vmask_e(1:jpi,1:jpj) = vmask(:,:,1) 
    241255      CALL lbc_lnk_icb( 'icbini', tmask_e, 'T', +1._wp, 1, 1 ) 
    242       CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 ) 
    243       CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 ) 
    244       ! 
     256      CALL lbc_lnk_icb( 'icbini', umask_e, 'U', +1._wp, 1, 1 ) 
     257      CALL lbc_lnk_icb( 'icbini', vmask_e, 'V', +1._wp, 1, 1 ) 
     258 
     259      ! definition of extended lat/lon array needed by icb_bilin_h 
     260      rlon_e(:,:) = 0._wp     ;  rlon_e(1:jpi,1:jpj) = glamt(:,:)  
     261      rlat_e(:,:) = 0._wp     ;  rlat_e(1:jpi,1:jpj) = gphit(:,:) 
     262      CALL lbc_lnk_icb( 'icbini', rlon_e, 'T', +1._wp, 1, 1 ) 
     263      CALL lbc_lnk_icb( 'icbini', rlat_e, 'T', +1._wp, 1, 1 ) 
     264      ! 
     265      ! definnitionn of extennded ff_f array needed by icb_utl_interp 
     266      ff_e(:,:) = 0._wp       ;  ff_e(1:jpi,1:jpj) = ff_f(:,:) 
     267      CALL lbc_lnk_icb( 'icbini', ff_e, 'F', +1._wp, 1, 1 ) 
     268 
    245269      ! assign each new iceberg with a unique number constructed from the processor number 
    246270      ! and incremented by the total number of processors 
    … …  
    338362               localpt%xi = REAL( mig(ji), wp ) 
    339363               localpt%yj = REAL( mjg(jj), wp ) 
    340                localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T' ) 
    341                localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T' ) 
     364               CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon )    
    342365               localpt%mass      = rn_initial_mass     (iberg) 
    343366               localpt%thickness = rn_initial_thickness(iberg) 
    … …  
    350373               localpt%uvel = 0._wp 
    351374               localpt%vvel = 0._wp 
     375               localpt%kb   = 1 
    352376               CALL icb_utl_incr() 
    353377               localberg%number(:) = num_bergs(:) 
    … …  
    383407         &              rn_bits_erosion_fraction        , rn_sicn_shift       , ln_passive_mode      ,   & 
    384408         &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   & 
    385          &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      ,                          & 
    386          &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out 
     409         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   & 
     410         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out        ,   & 
     411         &              ln_icb_grd 
    387412      !!---------------------------------------------------------------------- 
    388413 
    … …  
    463488            &                    'bits_erosion_fraction = ', rn_bits_erosion_fraction 
    464489 
     490         WRITE(numout,*) '   Use icb module modification from Merino et al. (2016) : ln_M2016 = ', ln_M2016 
     491         WRITE(numout,*) '       ground icebergs if icb bottom lvl hit the oce bottom level : ln_icb_grd = ', ln_icb_grd 
     492 
    465493         WRITE(numout,*) '   Shift of sea-ice concentration in erosion flux modulation ',   & 
    466494            &                    '(0<sicn_shift<1)    rn_sicn_shift  = ', rn_sicn_shift 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icbstp.F90

    r11536 r14054  
    5252CONTAINS 
    5353 
    54    SUBROUTINE icb_stp( kt ) 
     54   SUBROUTINE icb_stp( kt, Kmm ) 
    5555      !!---------------------------------------------------------------------- 
    5656      !!                  ***  ROUTINE icb_stp  *** 
    … …  
    6161      !!---------------------------------------------------------------------- 
    6262      INTEGER, INTENT(in) ::   kt   ! time step index 
     63      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    6364      ! 
    6465      LOGICAL ::   ll_sample_traj, ll_budget, ll_verbose   ! local logical 
    … …  
    7071      ! 
    7172      nktberg = kt 
     73      ! 
     74      !CALL test_icb_utl_getkb 
     75      !CALL ctl_stop('end test icb') 
    7276      ! 
    7377      IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data 
    … …  
    9296      ! 
    9397      !                                   !* copy nemo forcing arrays into iceberg versions with extra halo 
    94       CALL icb_utl_copy()                 ! only necessary for variables not on T points 
     98      CALL icb_utl_copy( Kmm )                 ! only necessary for variables not on T points 
    9599      ! 
    96100      ! 

  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/ICB/icbthm.F90

    r13281 r14054  
    4949      INTEGER, INTENT(in) ::   kt   ! timestep number, just passed to icb_utl_print_berg 
    5050      ! 
    51       INTEGER  ::   ii, ij 
    52       REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 
     51      INTEGER  ::   ii, ij, jk, ikb 
     52      REAL(wp) ::   zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb