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 13710 for NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL – NEMO

Ignore:
Timestamp:
2020-11-02T10:56:42+01:00 (4 years ago)
Author:
emanuelaclementi
Message:

branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves: merge with trunk@13708, see #2155 and #2339

Location:
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves
Files:
2 deleted
12 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13559        sette 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/EXPREF/context_nemo.xml

    r12276 r13710  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/EXPREF/file_def_nemo-oce.xml

    r9572 r13710  
    1515     <field field_ref="soce" />  
    1616     <field field_ref="ssh"  /> 
    17      <field field_ref="salgrad"  /> 
    18      <field field_ref="ke_zint"  /> 
     17     <field field_ref="socegrad"  /> 
     18     <field field_ref="eken_int"  /> 
    1919     <field field_ref="relvor"  /> 
    2020     <field field_ref="potvor"  /> 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/EXPREF/namelist_cfg

    r12489 r13710  
    2020&namusr_def    !   User defined :   CANAL configuration: Flat bottom, beta-plane 
    2121!----------------------------------------------------------------------- 
    22    rn_domszx   =   3600.   !  x horizontal size         [km] 
    23    rn_domszy   =   1800.   !  y horizontal size         [km] 
    24    rn_domszz   =   5000.   !  z vertical size            [m] 
    25    rn_dx       =     30.   !  x horizontal resolution   [km] 
    26    rn_dy       =     30.   !  y horizontal resolution   [km] 
    27    rn_dz       =    500.   !  z vertical resolution      [m] 
     22   rn_domszx   =   2000.   !  x horizontal size         [km] 
     23   rn_domszy   =   1000.   !  y horizontal size         [km] 
     24   rn_domszz   =   1000.   !  z vertical size            [m] 
     25   rn_dx       =     10.   !  x horizontal resolution   [km] 
     26   rn_dy       =     10.   !  y horizontal resolution   [km] 
     27   rn_dz       =   1000.   !  z vertical resolution      [m] 
    2828   rn_0xratio  =      0.5  !  x-domain ratio of the 0 
    2929   rn_0yratio  =      0.5  !  y-domain ratio of the 0 
    … …  
    3131   rn_ppgphi0  =    38.5   !  Reference latitude      [degrees] 
    3232   rn_u10      =      0.   !  10m wind speed              [m/s] 
    33      rn_windszx =   4000.   !  longitudinal wind extension   [km] 
    34      rn_windszy =   4000.   !  latitudinal wind extension    [km] 
    35      rn_uofac  =      0.   !  Uoce multiplicative factor (0.:absolute or 1.:relative winds) 
     33     rn_windszx =   90.    !  longitudinal wind extension   [km] 
     34     rn_windszy =   90.    !  latitudinal wind extension    [km] 
     35!!clem     rn_uofac  =     0.    !  Uoce multiplicative factor (0.:absolute or 1.:relative winds) 
    3636   rn_vtxmax   =      1.   !  initial vortex max current  [m/s] 
    3737   rn_uzonal   =      1.   !  initial zonal current       [m/s] 
    38      rn_ujetszx =   4000.   !  longitudinal jet extension   [km] 
    39      rn_ujetszy =   4000.   !  latitudinal jet extension    [km] 
     38     rn_ujetszx =   4000.  !  longitudinal jet extension   [km] 
     39     rn_ujetszy =   400.   !  latitudinal jet extension    [km] 
    4040   nn_botcase  =      0    !  bottom definition (0:flat, 1:bump) 
    41    nn_initcase =      1    !  initial condition case (0:rest, 1:zonal current, 2:current shear, 3: gaussian zonal current, 
    42       !                    !                          4: geostrophic zonal pulse, 5: vortex) 
    43    ln_sshnoise =  .false.  !  add random noise on initial ssh 
    44    rn_lambda   =     50.   ! gaussian lambda 
     41   nn_initcase =      1    !  initial condition case 
     42   !                       !          -1 : stratif at rest 
     43   !                       !           0 : rest 
     44   !                       !           1 : zonal current 
     45   !                       !           2 : current shear 
     46   !                       !           3 : gaussian zonal current 
     47   !                       !           4 : geostrophic zonal pulse 
     48   !                       !           5 : baroclinic vortex 
     49   ln_sshnoise =  .FALSE.  !  add random noise on initial ssh 
     50   rn_lambda   =     50.   !  gaussian lambda 
     51   nn_perio    = 1 
    4552/ 
    4653!----------------------------------------------------------------------- 
    … …  
    5966!----------------------------------------------------------------------- 
    6067   ln_linssh   =  .false.  !  =T  linear free surface  ==>>  model level are fixed in time 
    61    rn_Dt      =   1440.   !  time step for the dynamics (and tracer if nn_acc=0) 
    62    rn_atfp     =   0.05    !  asselin time filter parameter 
     68   rn_Dt      =   1200.    !  time step for the dynamics (and tracer if nn_acc=0) 
     69   rn_atfp     =   0.0     !  asselin time filter parameter 
     70/ 
     71!----------------------------------------------------------------------- 
     72&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg) 
     73!----------------------------------------------------------------------- 
     74   ln_write_cfg = .false.   !  (=T) create the domain configuration file 
     75      cn_domcfg_out = "domain_cfg" ! newly created domain configuration filename 
    6376/ 
    6477!!====================================================================== 
    … …  
    108121!!                                                                    !! 
    109122!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    110 !!   namdrg_top    top    friction                                      (ln_OFF =F & ln_isfcav=T) 
    111 !!   namdrg_bot    bottom friction                                      (ln_OFF =F) 
     123!!   namdrg_top    top    friction                                      (ln_drg_OFF =F & ln_isfcav=T) 
     124!!   namdrg_bot    bottom friction                                      (ln_drg_OFF =F) 
    112125!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    113126!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
    … …  
    117130&namdrg        !   top/bottom drag coefficient                          (default: NO selection) 
    118131!----------------------------------------------------------------------- 
    119    ln_OFF     = .true.    !  free-slip       : Cd = 0 
     132   ln_drg_OFF = .true.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
    120133/ 
    121134!!====================================================================== 
    … …  
    134147!----------------------------------------------------------------------- 
    135148   ln_seos     = .true.         !  = Use simplified equation of state (S-EOS) 
    136    !                            !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
     149   !                            !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
    137150   rn_a0       =  0.28        !  thermal expension coefficient (for simplified equation of state) 
    138151   rn_b0       =  0.          !  saline  expension coefficient (for simplified equation of state) 
    … …  
    148161   ln_traadv_OFF = .false. !  No tracer advection 
    149162   ln_traadv_cen = .false. !  2nd order centered scheme 
    150       nn_cen_h   =  4            !  =2/4, horizontal 2nd order CEN / 4th order CEN 
    151       nn_cen_v   =  4            !  =2/4, vertical   2nd order CEN / 4th order COMPACT 
     163      nn_cen_h   =  2            !  =2/4, horizontal 2nd order CEN / 4th order CEN 
     164      nn_cen_v   =  2            !  =2/4, vertical   2nd order CEN / 4th order COMPACT 
    152165   ln_traadv_fct = .false. !  FCT scheme 
    153       nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order 
     166      nn_fct_h   =  4            !  =2/4, horizontal 2nd / 4th order 
    154167      nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order 
    155168   ln_traadv_mus = .false. !  MUSCL scheme 
    … …  
    162175&namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection) 
    163176!----------------------------------------------------------------------- 
    164    ln_traldf_OFF   =  .true.  !  No explicit diffusion 
     177   !                       !  Operator type: 
     178   ln_traldf_OFF   = .true.    !  No explicit diffusion 
     179   ln_traldf_lap   = .false.   !    laplacian operator 
     180   ln_traldf_blp   = .false.   !  bilaplacian operator 
     181   ! 
     182   !                       !  Direction of action: 
     183   ln_traldf_lev   = .false.   !  iso-level 
     184   ln_traldf_hor   = .true.    !  horizontal  (geopotential) 
     185   ln_traldf_iso   = .false.   !  iso-neutral (standard operator) 
     186   ln_traldf_triad = .false.   !  iso-neutral (triad    operator) 
     187   ! 
     188   !                             !  iso-neutral options: 
     189   ln_traldf_msc   = .false.   !  Method of Stabilizing Correction      (both operators) 
     190   rn_slpmax       =  0.01     !  slope limit                           (both operators) 
     191   ln_triad_iso    = .false.   !  pure horizontal mixing in ML              (triad only) 
     192   rn_sw_triad     = 1         !  =1 switching triad ; =0 all 4 triads used (triad only) 
     193   ln_botmix_triad = .false.   !  lateral mixing on bottom                  (triad only) 
     194   ! 
     195   !                       !  Coefficients: 
     196   nn_aht_ijk_t    = 31         !  space/time variation of eddy coefficient: 
     197      !                             !   =-20 (=-30)    read in eddy_diffusivity_2D.nc (..._3D.nc) file 
     198      !                             !   =  0           constant 
     199      !                             !   = 10 F(k)      =ldf_c1d 
     200      !                             !   = 20 F(i,j)    =ldf_c2d 
     201      !                             !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation 
     202      !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d 
     203      !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
     204      !                        !  time invariant coefficients:  aht0 = 1/2  Ud*Ld   (lap case) 
     205      !                             !                           or   = 1/12 Ud*Ld^3 (blp case) 
     206      rn_Ud        = 0.01           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) 
     207      rn_Ld        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10) 
    165208/ 
    166209!!====================================================================== 
    … …  
    183226      nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction 
    184227   ln_dynadv_cen2 = .false. !  flux form - 2nd order centered scheme 
    185    ln_dynadv_ubs = .true.  !  flux form - 3rd order UBS      scheme 
     228   ln_dynadv_ubs  = .true.  !  flux form - 3rd order UBS      scheme 
    186229/ 
    187230!----------------------------------------------------------------------- 
    188231&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    189232!----------------------------------------------------------------------- 
    190    ln_dynvor_ene = .true.  !  energy conserving scheme 
    191    ln_dynvor_ens = .false. !  enstrophy conserving scheme 
    192    ln_dynvor_mix = .false. !  mixed scheme 
     233   ln_dynvor_ene = .false.  !  energy conserving scheme 
     234   ln_dynvor_ens = .false.  !  enstrophy conserving scheme 
     235   ln_dynvor_mix = .false.  !  mixed scheme 
    193236   ln_dynvor_een = .false.  !  energy & enstrophy scheme 
     237   ln_dynvor_enT = .false.  !  energy conserving scheme (T-point) 
     238   ln_dynvor_eeT = .true.   !  energy conserving scheme (een using e3t) 
    194239      nn_een_e3f = 0             !  e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    195240/ 
    … …  
    210255         !                          !                     = 1 Boxcar over   nn_e sub-steps 
    211256         !                          !                     = 2 Boxcar over 2*nn_e  "    " 
    212       ln_bt_auto    = .false.    ! Number of sub-step defined from: 
     257      ln_bt_auto    = .true.    ! Number of sub-step defined from: 
    213258         nn_e      =  24         ! =F : the number of sub-step in rn_Dt seconds 
    214259/ 
    … …  
    222267   !                       !  Direction of action  : 
    223268   ln_dynldf_lev =  .false.    !  iso-level 
    224    ln_dynldf_hor =  .true.    !  horizontal (geopotential) 
     269   ln_dynldf_hor =  .false.    !  horizontal (geopotential) 
    225270   ln_dynldf_iso =  .false.    !  iso-neutral 
    226271   !                       !  Coefficient 
    227    nn_ahm_ijk_t  = 20           !  space/time variation of eddy coef 
     272   nn_ahm_ijk_t  = 31           !  space/time variation of eddy coef 
    228273      !                             !  =-30  read in eddy_viscosity_3D.nc file 
    229274      !                             !  =-20  read in eddy_viscosity_2D.nc file 
    … …  
    271316!!                                                                    !! 
    272317!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    273 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    274318!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    275319!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    276320!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
     321<<<<<<< .working 
    277322!!   namflo       float parameters                                      (default: OFF) 
    278323!!   nam_diadct   transports through some sections                      (default: OFF) 
     324||||||| .merge-left.r13465 
     325!!   namflo       float parameters                                      (default: OFF) 
     326!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF) 
     327!!   nam_diadct   transports through some sections                      (default: OFF) 
     328======= 
     329!!   namflo       float parameters                                      ("key_float") 
     330!!   nam_diaharm  Harmonic analysis of tidal constituents               ("key_diaharm") 
     331!!   namdct       transports through some sections                      ("key_diadct") 
     332!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF) 
     333>>>>>>> .merge-right.r13470 
    279334!!   nam_dia25h   25h Mean Output                                       (default: OFF) 
    280335!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4") 
    … …  
    285340!----------------------------------------------------------------------- 
    286341   ln_glo_trd  = .false.   ! (T) global domain averaged diag for T, T^2, KE, and PE 
    287    ln_dyn_trd  = .true.   ! (T) 3D momentum trend output 
     342   ln_dyn_trd  = .true.    ! (T) 3D momentum trend output 
    288343   ln_dyn_mxl  = .false.   ! (T) 2D momentum trends averaged over the mixed layer (not coded yet) 
    289344   ln_vor_trd  = .false.   ! (T) 2D barotropic vorticity trends (not coded yet) 
    … …  
    312367&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi") 
    313368!----------------------------------------------------------------------- 
     369!!   jpni        =   8       !  jpni   number of processors following i (set automatically if < 1) 
     370!!   jpnj        =   1       !  jpnj   number of processors following j (set automatically if < 1) 
    314371/ 
    315372!----------------------------------------------------------------------- 
    316373&namctl        !   Control prints                                       (default: OFF) 
    317374!----------------------------------------------------------------------- 
     375   ln_timing   = .true.   !  timing by routine write out in timing.output file 
     376!!   ln_diacfl   = .true.   !  CFL diagnostics write out in cfl_diagnostics.ascii 
    318377/ 
    319378!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/domvvl.F90

    r12489 r13710  
    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 
    1112   !!---------------------------------------------------------------------- 
    1213 
    13    !!---------------------------------------------------------------------- 
    14    !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
    15    !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
    16    !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
    17    !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
    18    !!   dom_vvl_rst      : read/write restart file 
    19    !!   dom_vvl_ctl      : Check the vvl options 
    20    !!---------------------------------------------------------------------- 
    2114   USE oce             ! ocean dynamics and tracers 
    2215   USE phycst          ! physical constant 
    … …  
    3629   PRIVATE 
    3730 
    38    PUBLIC  dom_vvl_init       ! called by domain.F90 
    39    PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    40    PUBLIC  dom_vvl_sf_update  ! called by step.F90 
    41    PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    42  
    4331   !                                                      !!* Namelist nam_vvl 
    4432   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate 
    … …  
    6250   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6351 
     52#if defined key_qco 
     53   !!---------------------------------------------------------------------- 
     54   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate 
     55   !!---------------------------------------------------------------------- 
     56#else 
     57   !!---------------------------------------------------------------------- 
     58   !!   Default key      Old management of time varying vertical coordinate 
     59   !!---------------------------------------------------------------------- 
     60    
     61   !!---------------------------------------------------------------------- 
     62   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
     63   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
     64   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
     65   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
     66   !!   dom_vvl_rst      : read/write restart file 
     67   !!   dom_vvl_ctl      : Check the vvl options 
     68   !!---------------------------------------------------------------------- 
     69 
     70   PUBLIC  dom_vvl_init       ! called by domain.F90 
     71   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
     72   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
     73   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
     74   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
     75    
     76   !! * Substitutions 
     77#  include "do_loop_substitute.h90" 
    6478   !!---------------------------------------------------------------------- 
    6579   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    116130      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
    117131      ! 
     132      IF(lwp) WRITE(numout,*) 
     133      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
     134      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     135      ! 
     136      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
     137      ! 
     138      !                    ! Allocate module arrays 
     139      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
     140      ! 
     141      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     142      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     143      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     144      ! 
     145      CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
     146      ! 
     147   END SUBROUTINE dom_vvl_init 
     148 
     149 
     150   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 
     151      !!---------------------------------------------------------------------- 
     152      !!                ***  ROUTINE dom_vvl_init  *** 
     153      !!                    
     154      !! ** Purpose :  Interpolation of all scale factors,  
     155      !!               depths and water column heights 
     156      !! 
     157      !! ** Method  :  - interpolate scale factors 
     158      !! 
     159      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     160      !!              - Regrid: e3(u/v)_n 
     161      !!                        e3(u/v)_b        
     162      !!                        e3w_n            
     163      !!                        e3(u/v)w_b       
     164      !!                        e3(u/v)w_n       
     165      !!                        gdept_n, gdepw_n and gde3w_n 
     166      !!              - h(t/u/v)_0 
     167      !!              - frq_rst_e3t and frq_rst_hdv 
     168      !! 
     169      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     170      !!---------------------------------------------------------------------- 
     171      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     172      !!---------------------------------------------------------------------- 
    118173      INTEGER ::   ji, jj, jk 
    119174      INTEGER ::   ii0, ii1, ij0, ij1 
    120175      REAL(wp)::   zcoef 
    121176      !!---------------------------------------------------------------------- 
    122       ! 
    123       IF(lwp) WRITE(numout,*) 
    124       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
    125       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    126       ! 
    127       CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    128       ! 
    129       !                    ! Allocate module arrays 
    130       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    131       ! 
    132       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    133       CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
    134       e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    135177      ! 
    136178      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    … …  
    160202      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
    161203      gdepw(:,:,1,Kbb) = 0.0_wp 
    162       DO jk = 2, jpk                               ! vertical sum 
    163          DO jj = 1,jpj 
    164             DO ji = 1,jpi 
    165                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    166                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    167                !                             ! 0.5 where jk = mikt      
     204      DO_3D( 1, 1, 1, 1, 2, jpk ) 
     205         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     206         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     207         !                             ! 0.5 where jk = mikt      
    168208!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    169                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    170                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    171                gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    172                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
    173                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    174                gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    175                gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    176                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    177             END DO 
    178          END DO 
    179       END DO 
     209         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
     210         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     211         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     212            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     213         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     214         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     215         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     216            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     217      END_3D 
    180218      ! 
    181219      !                    !==  thickness of the water column  !!   (ocean portion only) 
    … …  
    212250         ENDIF 
    213251         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    214             DO jj = 1, jpj 
    215                DO ji = 1, jpi 
     252            DO_2D( 1, 1, 1, 1 ) 
    216253!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    217                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    218                      ! values outside the equatorial band and transition zone (ztilde) 
    219                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    220                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    221                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    222                      ! values inside the equatorial band (ztilde as zstar) 
    223                      frq_rst_e3t(ji,jj) =  0.0_wp 
    224                      frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
    225                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    226                      !                                      ! (linearly transition from z-tilde to z-star) 
    227                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    228                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    229                         &                                          * 180._wp / 3.5_wp ) ) 
    230                      frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
    231                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
    232                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    233                         &                                          * 180._wp / 3.5_wp ) ) 
    234                   ENDIF 
    235                END DO 
    236             END DO 
     254               IF( ABS(gphit(ji,jj)) >= 6.) THEN 
     255                  ! values outside the equatorial band and transition zone (ztilde) 
     256                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
     257                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     258               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     259                  ! values inside the equatorial band (ztilde as zstar) 
     260                  frq_rst_e3t(ji,jj) =  0.0_wp 
     261                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
     262               ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
     263                  !                                      ! (linearly transition from z-tilde to z-star) 
     264                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
     265                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     266                     &                                          * 180._wp / 3.5_wp ) ) 
     267                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
     268                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
     269                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     270                     &                                          * 180._wp / 3.5_wp ) ) 
     271               ENDIF 
     272            END_2D 
    237273            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    238274               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
    239                   ii0 = 103   ;   ii1 = 111        
    240                   ij0 = 128   ;   ij1 = 135   ;    
     275                  ii0 = 103 + nn_hls - 1   ;   ii1 = 111 + nn_hls - 1       
     276                  ij0 = 128 + nn_hls       ;   ij1 = 135 + nn_hls 
    241277                  frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
    242278                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rn_Dt 
    … …  
    264300      ENDIF 
    265301      ! 
    266    END SUBROUTINE dom_vvl_init 
     302   END SUBROUTINE dom_vvl_zgr 
    267303 
    268304 
    … …  
    298334      LOGICAL                ::   ll_do_bclinic         ! local logical 
    299335      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv 
    300       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t 
     336      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3t 
     337      LOGICAL , DIMENSION(:,:,:), ALLOCATABLE ::   llmsk 
    301338      !!---------------------------------------------------------------------- 
    302339      ! 
    … …  
    329366      END DO 
    330367      ! 
    331       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    332          !                                                            ! ------baroclinic part------ ! 
     368      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     369         !                                                               ! ------baroclinic part------ ! 
    333370         ! I - initialization 
    334371         ! ================== 
    … …  
    383420         zwu(:,:) = 0._wp 
    384421         zwv(:,:) = 0._wp 
    385          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    386             DO jj = 1, jpjm1 
    387                DO ji = 1, fs_jpim1   ! vector opt. 
    388                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    389                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    390                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    391                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    392                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    393                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    394                END DO 
    395             END DO 
    396          END DO 
    397          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    398             DO ji = 1, jpi 
    399                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    400                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    401             END DO 
    402          END DO 
    403          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    404             DO jj = 2, jpjm1 
    405                DO ji = fs_2, fs_jpim1   ! vector opt. 
    406                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    407                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    408                      &                                            ) * r1_e1e2t(ji,jj) 
    409                END DO 
    410             END DO 
    411          END DO 
     422         DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     423            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
     424               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     425            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     426               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     427            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     428            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     429         END_3D 
     430         DO_2D( 1, 1, 1, 1 ) 
     431            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     432            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     433         END_2D 
     434         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     435            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     436               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     437               &                                            ) * r1_e1e2t(ji,jj) 
     438         END_3D 
    412439         !                       ! d - thickness diffusion transport: boundary conditions 
    413440         !                             (stored for tracer advction and continuity equation) 
    … …  
    416443         ! 4 - Time stepping of baroclinic scale factors 
    417444         ! --------------------------------------------- 
    418          ! Leapfrog time stepping 
    419          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    420445         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    421446         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
    … …  
    423448         ! Maximum deformation control 
    424449         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    425          ze3t(:,:,jpk) = 0._wp 
    426          DO jk = 1, jpkm1 
    427             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
    428          END DO 
    429          z_tmax = MAXVAL( ze3t(:,:,:) ) 
    430          CALL mpp_max( 'domvvl', z_tmax )                 ! max over the global domain 
    431          z_tmin = MINVAL( ze3t(:,:,:) ) 
    432          CALL mpp_min( 'domvvl', z_tmin )                 ! min over the global domain 
     450         ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) ) 
     451         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     452            ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     453         END_3D 
     454         ! 
     455         llmsk(   1:Nis1,:,:) = .FALSE.   ! exclude halos from the checked region 
     456         llmsk(Nie1: jpi,:,:) = .FALSE. 
     457         llmsk(:,   1:Njs1,:) = .FALSE. 
     458         llmsk(:,Nje1: jpj,:) = .FALSE. 
     459         ! 
     460         llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp                  ! define only the inner domain 
     461         z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_max( 'domvvl', z_tmax )   ! max over the global domain 
     462         z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_min( 'domvvl', z_tmin )   ! min over the global domain 
    433463         ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    434464         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
    435             IF( lk_mpp ) THEN 
    436                CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max ) 
    437                CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min ) 
    438             ELSE 
    439                ijk_max = MAXLOC( ze3t(:,:,:) ) 
    440                ijk_max(1) = ijk_max(1) + nimpp - 1 
    441                ijk_max(2) = ijk_max(2) + njmpp - 1 
    442                ijk_min = MINLOC( ze3t(:,:,:) ) 
    443                ijk_min(1) = ijk_min(1) + nimpp - 1 
    444                ijk_min(2) = ijk_min(2) + njmpp - 1 
    445             ENDIF 
     465            CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max ) 
     466            CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min ) 
    446467            IF (lwp) THEN 
    447468               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
    … …  
    452473            ENDIF 
    453474         ENDIF 
     475         DEALLOCATE( ze3t, llmsk ) 
    454476         ! - ML - end test 
    455477         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    … …  
    613635         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    614636      ENDIF 
    615       gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
    616       gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    617  
    618       e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 
    619       e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 
    620       e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 
    621637 
    622638      ! Compute all missing vertical scale factor and depths 
    … …  
    641657      gdepw(:,:,1,Kmm) = 0.0_wp 
    642658      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    643       DO jk = 2, jpk 
    644          DO jj = 1,jpj 
    645             DO ji = 1,jpi 
    646               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    647                                                                  ! 1 for jk = mikt 
    648                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    649                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    650                gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    651                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
    652                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    653             END DO 
    654          END DO 
    655       END DO 
     659      DO_3D( 1, 1, 1, 1, 2, jpk ) 
     660        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     661                                                           ! 1 for jk = mikt 
     662         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     663         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     664         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     665             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     666         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     667      END_3D 
    656668 
    657669      ! Local depth and Inverse of the local depth of the water 
    … …  
    700712         ! 
    701713      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    702          DO jk = 1, jpk 
    703             DO jj = 1, jpjm1 
    704                DO ji = 1, fs_jpim1   ! vector opt. 
    705                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    706                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    707                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    708                END DO 
    709             END DO 
    710          END DO 
     714         DO_3D( 1, 0, 1, 0, 1, jpk ) 
     715            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     716               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     717               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     718         END_3D 
    711719         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    712720         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    713721         ! 
    714722      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    715          DO jk = 1, jpk 
    716             DO jj = 1, jpjm1 
    717                DO ji = 1, fs_jpim1   ! vector opt. 
    718                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    719                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    720                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    721                END DO 
    722             END DO 
    723          END DO 
     723         DO_3D( 1, 0, 1, 0, 1, jpk ) 
     724            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     725               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     726               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     727         END_3D 
    724728         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    725729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    726730         ! 
    727731      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    728          DO jk = 1, jpk 
    729             DO jj = 1, jpjm1 
    730                DO ji = 1, fs_jpim1   ! vector opt. 
    731                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    732                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    733                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    734                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    735                END DO 
    736             END DO 
    737          END DO 
     732         DO_3D( 1, 0, 1, 0, 1, jpk ) 
     733            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     734               &                       *    r1_e1e2f(ji,jj)                                                  & 
     735               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     736               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     737         END_3D 
    738738         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    739739         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
    … …  
    803803         IF( ln_rstart ) THEN                   !* Read the restart file 
    804804            CALL rst_read_open                  !  open the restart file if necessary 
    805             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
     805            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
    806806            ! 
    807807            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
    … …  
    810810            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    811811            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     812            ! 
    812813            !                             ! --------- ! 
    813814            !                             ! all cases ! 
    814815            !                             ! --------- ! 
     816            ! 
    815817            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    816                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    817                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
     818               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     819               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    818820               ! needed to restart if land processor not computed  
    819821               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    … …  
    828830               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    829831               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    830                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
    831                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     832               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
     833               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    832834               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    833835               l_1st_euler = .true. 
    … …  
    835837               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    836838               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    837                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
    838                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
     839               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
     840               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    839841               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    840842               l_1st_euler = .true. 
    … …  
    842844               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    843845               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    844                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     846               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    845847               DO jk = 1, jpk 
    846848                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
    … …  
    861863               !                          ! ----------------------- ! 
    862864               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
    863                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios ) 
    864                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios ) 
     865                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios ) 
     866                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios ) 
    865867               ELSE                            ! one at least array is missing 
    866868                  tilde_e3t_b(:,:,:) = 0.0_wp 
    … …  
    871873                  !                       ! ------------ ! 
    872874                  IF( id5 > 0 ) THEN  ! required array exists 
    873                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 
     875                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 
    874876                  ELSE                ! array is missing 
    875877                     hdiv_lf(:,:,:) = 0.0_wp 
    … …  
    895897                  ssh(:,:,Kbb) = -ssh_ref 
    896898 
    897                   DO jj = 1, jpj 
    898                      DO ji = 1, jpi 
    899                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    900                            ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    901                            ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    902                         ENDIF 
    903                      ENDDO 
    904                   ENDDO 
     899                  DO_2D( 1, 1, 1, 1 ) 
     900                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     901                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     902                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     903                     ENDIF 
     904                  END_2D 
    905905               ENDIF !If test case else 
    906906 
    … …  
    913913               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    914914 
    915                DO ji = 1, jpi 
    916                   DO jj = 1, jpj 
    917                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    918                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    919                      ENDIF 
    920                   END DO  
    921                END DO  
     915               DO_2D( 1, 1, 1, 1 ) 
     916                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     917                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     918                  ENDIF 
     919               END_2D 
    922920               ! 
    923921            ELSE 
    924922               ! 
    925                ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 
    926                CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    927                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 
     923               ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 
     924               ! 
     925               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )   
     926               ! 
     927               ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 
    928928               ! 
    929929               DO jk=1,jpk 
    930                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    931                     &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    932                     &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points 
     930                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
     931                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     932                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points 
    933933               END DO 
    934934               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    935                ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w 
     935               ssh(:,:,Kmm) = ssh(:,:,Kbb)                                     ! needed later for gde3w 
    936936               ! 
    937937            END IF           ! end of ll_wd edits 
    … …  
    10251025      ! 
    10261026      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    1027       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    10281027      ! 
    10291028      IF(lwp) THEN                   ! Print the choice 
    … …  
    10411040   END SUBROUTINE dom_vvl_ctl 
    10421041 
     1042#endif 
     1043 
    10431044   !!====================================================================== 
    10441045END MODULE domvvl 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/stpctl.F90

    r12377 r13710  
    1919   USE dom_oce         ! ocean space and time domain variables  
    2020   USE c1d             ! 1D vertical configuration 
     21   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
     22   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
     23   !   
    2124   USE diawri          ! Standard run outputs       (dia_wri_state routine) 
    22    ! 
    2325   USE in_out_manager  ! I/O manager 
    2426   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2527   USE lib_mpp         ! distributed memory computing 
    26    USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
    27    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
    28  
     28   ! 
    2929   USE netcdf          ! NetCDF library 
    3030   IMPLICIT NONE 
    … …  
    3333   PUBLIC stp_ctl           ! routine called by step.F90 
    3434 
    35    INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 
    36    LOGICAL  ::   lsomeoce 
     35   INTEGER                ::   nrunid   ! netcdf file id 
     36   INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id 
    3737   !!---------------------------------------------------------------------- 
    3838   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4242CONTAINS 
    4343 
    44    SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 
     44   SUBROUTINE stp_ctl( kt, Kmm ) 
    4545      !!---------------------------------------------------------------------- 
    4646      !!                    ***  ROUTINE stp_ctl  *** 
    … …  
    4949      !! 
    5050      !! ** Method  : - Save the time step in numstp 
    51       !!              - Print it each 50 time steps 
    52       !!              - Stop the run IF problem encountered by setting indic=-3 
     51      !!              - Stop the run IF problem encountered by setting nstop > 0 
    5352      !!                Problems checked: |ssh| maximum larger than 10 m 
    5453      !!                                  |U|   maximum larger than 10 m/s  
    … …  
    5756      !! ** Actions :   "time.step" file = last ocean time-step 
    5857      !!                "run.stat"  file = run statistics 
    59       !!                nstop indicator sheared among all local domain (lk_mpp=T) 
     58      !!                 nstop indicator sheared among all local domain 
    6059      !!---------------------------------------------------------------------- 
    6160      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    62       INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index 
    63       INTEGER, INTENT(inout) ::   kindic   ! error indicator 
    64       !! 
    65       INTEGER                ::   ji, jj, jk          ! dummy loop indices 
    66       INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices 
    67       INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
    68       REAL(wp)               ::   zzz                 ! local real  
    69       REAL(wp), DIMENSION(9) ::   zmax 
    70       LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
    71       CHARACTER(len=20) :: clname 
    72       !!---------------------------------------------------------------------- 
    73       ! 
    74       ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
    75       ll_colruns = ll_wrtstp .AND. ( ln_ctl .OR. sn_cfctl%l_runstat ) 
    76       ll_wrtruns = ll_colruns .AND. lwm 
    77       IF( kt == nit000 .AND. lwp ) THEN 
    78          WRITE(numout,*) 
    79          WRITE(numout,*) 'stp_ctl : time-stepping control' 
    80          WRITE(numout,*) '~~~~~~~' 
    81          !                                ! open time.step file 
    82          IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    83          !                                ! open run.stat file(s) at start whatever 
    84          !                                ! the value of sn_cfctl%ptimincr 
    85          IF( lwm .AND. ( ln_ctl .OR. sn_cfctl%l_runstat ) ) THEN 
     61      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index 
     62      !! 
     63      INTEGER                         ::   ji                                    ! dummy loop indices 
     64      INTEGER                         ::   idtime, istatus 
     65      INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax 
     66      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices 
     67      REAL(wp)                        ::   zzz                                   ! local real  
     68      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
     69      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 
     70      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
     71      CHARACTER(len=20)               ::   clname 
     72      !!---------------------------------------------------------------------- 
     73      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid 
     74      ! 
     75      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
     76      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1  
     77      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 
     78      ! 
     79      IF( kt == nit000 ) THEN 
     80         ! 
     81         IF( lwp ) THEN 
     82            WRITE(numout,*) 
     83            WRITE(numout,*) 'stp_ctl : time-stepping control' 
     84            WRITE(numout,*) '~~~~~~~' 
     85         ENDIF 
     86         !                                ! open time.step    ascii file, done only by 1st subdomain 
     87         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     88         ! 
     89         IF( ll_wrtruns ) THEN 
     90            !                             ! open run.stat     ascii file, done only by 1st subdomain 
    8691            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     92            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain 
    8793            clname = 'run.stat.nc' 
    8894            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
    89             istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun ) 
    90             istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime ) 
    91             istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 
    92             istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu  ) 
    93             istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1 ) 
    94             istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2 ) 
    95             istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1 ) 
    96             istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2 ) 
     95            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid ) 
     96            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime ) 
     97            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) ) 
     98            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) ) 
     99            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) ) 
     100            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) ) 
     101            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) ) 
     102            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) ) 
    97103            IF( ln_zad_Aimp ) THEN 
    98                istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 
    99                istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1 ) 
     104               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) ) 
     105               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) ) 
    100106            ENDIF 
    101             istatus = NF90_ENDDEF(idrun) 
    102             zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use 
    103          ENDIF 
    104       ENDIF 
    105       IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 
    106       ! 
    107       IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file) 
     107            istatus = NF90_ENDDEF(nrunid) 
     108         ENDIF 
     109         !     
     110      ENDIF 
     111      ! 
     112      !                                   !==              write current time step              ==! 
     113      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
     114      IF( lwm .AND. ll_wrtstp ) THEN 
    108115         WRITE ( numstp, '(1x, i8)' )   kt 
    109116         REWIND( numstp ) 
    110117      ENDIF 
    111       ! 
    112       !                                   !==  test of extrema  ==! 
     118      !                                   !==            test of local extrema           ==! 
     119      !                                   !==  done by all processes at every time step  ==! 
     120      ! 
     121      llmsk(   1:Nis1,:,:) = .FALSE.                                              ! exclude halos from the checked region 
     122      llmsk(Nie1: jpi,:,:) = .FALSE. 
     123      llmsk(:,   1:Njs1,:) = .FALSE. 
     124      llmsk(:,Nje1: jpj,:) = .FALSE. 
     125      ! 
     126      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain 
     127      ! 
     128      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain? 
     129      ! 
    113130      IF( ll_wd ) THEN 
    114          zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
     131         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
    115132      ELSE 
    116          zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max 
    117       ENDIF 
    118       zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
    119       zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
    120       zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
    121       zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
    122       zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
    123       zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
    124       IF( ln_zad_Aimp ) THEN 
    125          zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
    126          zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max 
    127       ENDIF 
    128       ! 
     133         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max 
     134      ENDIF 
     135      llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     136      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only) 
     137      llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     138      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max 
     139      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max 
     140      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file 
     141         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max 
     142         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max 
     143         IF( ln_zad_Aimp ) THEN 
     144            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max 
     145            llmsk(:,:,:) = wmask(:,:,:) == 1._wp 
     146            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max 
     147         ELSE 
     148            zmax(7:8) = 0._wp 
     149         ENDIF 
     150      ELSE 
     151         zmax(5:8) = 0._wp 
     152      ENDIF 
     153      zmax(9) = REAL( nstop, wp )                                                 ! stop indicator 
     154      ! 
     155      !                                   !==               get global extrema             ==! 
     156      !                                   !==  done by all processes if writting run.stat  ==! 
    129157      IF( ll_colruns ) THEN 
    130          CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    131          nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains 
    132       ENDIF 
    133       !                                   !==  run statistics  ==!   ("run.stat" files) 
     158         zmaxlocal(:) = zmax(:) 
     159         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true.  
     160         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
     161      ELSE 
     162         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 
     163         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values... 
     164      ENDIF 
     165      ! 
     166      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage!  
     167      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage! 
     168      IF( ll_colruns ) THEN 
     169         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage!  
     170         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage! 
     171      ENDIF 
     172      ! 
     173      !                                   !==              write "run.stat" files              ==! 
     174      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    134175      IF( ll_wrtruns ) THEN 
    135          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
    136          istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 
    137          istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) ) 
    138          istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) ) 
    139          istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) ) 
    140          istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) ) 
    141          istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) ) 
    142          IF( ln_zad_Aimp ) THEN 
    143             istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) ) 
    144             istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) ) 
    145          ENDIF 
    146          IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
    147          IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
     176         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 
     177         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 
     178            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 
     179         END DO 
     180         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
    148181      END IF 
    149       !                                   !==  error handling  ==! 
    150       IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 
    151          &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
    152          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    153 !!$      &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
    154 !!$      &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
    155 !!$      &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
    156          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
    157          IF( lk_mpp .AND. ln_ctl ) THEN 
    158             CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
    159             CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
    160             CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 
    161             CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 
     182      !                                   !==               error handling               ==! 
     183      !                                   !==  done by all processes at every time step  ==! 
     184      ! 
     185      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
     186         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
     187!!$         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
     188!!$         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
     189!!$         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
     190         & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
     191         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     192         ! 
     193         iloc(:,:) = 0 
     194         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 
     195            ! first: close the netcdf file, so we can read it 
     196            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid) 
     197            ! get global loc on the min/max 
     198            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain 
     199            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), llmsk(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F  
     200            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     201            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)), llmsk(:,:,:), zzz, iloc(1:3,2) ) 
     202            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     203            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,3) ) 
     204            CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,4) ) 
     205            ! find which subdomain has the max. 
     206            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0 
     207            DO ji = 1, 9 
     208               IF( zmaxlocal(ji) == zmax(ji) ) THEN 
     209                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1 
     210               ENDIF 
     211            END DO 
     212            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain 
     213            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain 
     214            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain 
     215         ELSE                    ! find local min and max locations: 
     216            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 
     217            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain 
     218            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) ) 
     219            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     220            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask = llmsk(:,:,:) ) 
     221            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     222            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) ) 
     223            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) ) 
     224            DO ji = 1, 4   ! local domain indices ==> global domain indices, excluding halos 
     225               iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /) 
     226            END DO 
     227            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information 
     228         ENDIF 
     229         ! 
     230         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
     231         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
     232         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
     233         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
     234         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
     235         IF( Agrif_Root() ) THEN 
     236            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 
    162237         ELSE 
    163             ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
    164             iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    165             is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    166             is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    167          ENDIF 
    168           
    169          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    170          WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
    171          WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
    172          WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 
    173          WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3) 
    174          WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file' 
    175           
     238            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 
     239         ENDIF 
     240         ! 
    176241         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
    177           
    178          IF( .NOT. ln_ctl ) THEN 
    179             WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 
    180             CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 ) 
    181          ELSE 
    182             CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' ) 
    183          ENDIF 
    184  
    185          kindic = -3 
    186          ! 
    187       ENDIF 
    188       ! 
    189 9100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5) 
    190 9200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
    191 9300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5) 
    192 9400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5) 
     242         ! 
     243         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files 
     244            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
     245            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop) 
     246            ENDIF 
     247         ELSE                                    ! only mpi subdomains with errors are here -> STOP now 
     248            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
     249         ENDIF 
     250         ! 
     251      ENDIF 
     252      ! 
     253      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet... 
     254         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error 
     255         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock 
     256      ENDIF 
     257      ! 
    1932589500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
    194259      ! 
    195260   END SUBROUTINE stp_ctl 
     261 
     262 
     263   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax ) 
     264      !!---------------------------------------------------------------------- 
     265      !!                     ***  ROUTINE wrt_line  *** 
     266      !! 
     267      !! ** Purpose :   write information line 
     268      !! 
     269      !!---------------------------------------------------------------------- 
     270      CHARACTER(len=*),      INTENT(  out) ::   cdline 
     271      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix 
     272      REAL(wp),              INTENT(in   ) ::   pval 
     273      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc 
     274      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax 
     275      ! 
     276      CHARACTER(len=80) ::   clsuff 
     277      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax 
     278      CHARACTER(len=9 ) ::   cli, clj, clk 
     279      CHARACTER(len=1 ) ::   clfmt 
     280      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why 
     281      INTEGER           ::   ifmtk 
     282      !!---------------------------------------------------------------------- 
     283      WRITE(clkt , '(i9)') kt 
     284       
     285      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9) 
     286      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF 
     287      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum 
     288      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9) 
     289      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1 
     290                                   WRITE(clmax, cl4) kmax-1 
     291      ! 
     292      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9) 
     293      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF 
     294      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9) 
     295      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF 
     296      ! 
     297      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin) 
     298      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax) 
     299      ENDIF 
     300      IF(kloc(3) == 0) THEN 
     301         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9) 
     302         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string 
     303         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff) 
     304      ELSE 
     305         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9) 
     306         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF 
     307         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF 
     308         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff) 
     309      ENDIF 
     310      ! 
     3119100  FORMAT('MPI rank ', a) 
     3129200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 
     3139300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a) 
     3149400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 
     315      ! 
     316   END SUBROUTINE wrt_line 
     317 
    196318 
    197319   !!====================================================================== 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/trazdf.F90

    r12489 r13710  
    3535   PUBLIC   tra_zdf_imp   ! called by trczdf.F90 
    3636 
     37   !! * Substitutions 
     38#  include "do_loop_substitute.h90" 
    3739   !!---------------------------------------------------------------------- 
    3840   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7779      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    7880      ! JMM : restore negative salinities to small salinities: 
    79 !!$   WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
     81!!$      WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
    8082!!gm 
    8183 
    … …  
    9597      ENDIF 
    9698      !                                          ! print mean trends (used for debugging) 
    97       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
    98          &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     99      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
     100         &                                  tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    99101      ! 
    100102      IF( ln_timing )   CALL timing_stop('tra_zdf') 
    … …  
    154156            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
    155157               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
    156                   DO jk = 2, jpkm1 
    157                      DO jj = 2, jpjm1 
    158                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    159                            zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
    160                         END DO 
    161                      END DO 
    162                   END DO 
     158                  DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     159                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     160                  END_3D 
    163161               ELSE                          ! standard or triad iso-neutral operator 
    164                   DO jk = 2, jpkm1 
    165                      DO jj = 2, jpjm1 
    166                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                            zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
    168                         END DO 
    169                      END DO 
    170                   END DO 
     162                  DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     163                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     164                  END_3D 
    171165               ENDIF 
    172166            ENDIF 
    … …  
    174168            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    175169            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection 
    176                DO jk = 1, jpkm1 
    177                   DO jj = 2, jpjm1 
    178                      DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
    179                         zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
    180                         zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
    181                         zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
    182                            &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
    183                         zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
    184                         zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp ) 
    185                     END DO 
    186                   END DO 
    187                END DO 
     170               DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     171                  zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
     172                  zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     173                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
     174                     &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
     175                  zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
     176                  zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp ) 
     177               END_3D 
    188178            ELSE 
    189                DO jk = 1, jpkm1 
    190                   DO jj = 2, jpjm1 
    191                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    192                         zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
    193                         zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
    194                         zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    195                     END DO 
    196                   END DO 
    197                END DO 
     179               DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     180                  zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
     181                  zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     182                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     183               END_3D 
    198184            ENDIF 
    199185            ! 
    … …  
    217203            !   used as a work space array: its value is modified. 
    218204            ! 
    219             DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    220                DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction) 
    221                   zwt(ji,jj,1) = zwd(ji,jj,1) 
    222                END DO 
    223             END DO 
    224             DO jk = 2, jpkm1 
    225                DO jj = 2, jpjm1 
    226                   DO ji = fs_2, fs_jpim1 
    227                      zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    228                   END DO 
    229                END DO 
    230             END DO 
     205            DO_2D( 0, 0, 0, 0 ) 
     206               zwt(ji,jj,1) = zwd(ji,jj,1) 
     207            END_2D 
     208            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     209               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     210            END_3D 
    231211            ! 
    232212         ENDIF  
    233213         !          
    234          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    235             DO ji = fs_2, fs_jpim1 
    236                pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
    237             END DO 
    238          END DO 
    239          DO jk = 2, jpkm1 
    240             DO jj = 2, jpjm1 
    241                DO ji = fs_2, fs_jpim1 
    242                   zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
    243                   pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
    244                END DO 
    245             END DO 
    246          END DO 
     214         DO_2D( 0, 0, 0, 0 ) 
     215            pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
     216         END_2D 
     217         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     218            zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
     219            pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
     220         END_3D 
    247221         ! 
    248          DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    249             DO ji = fs_2, fs_jpim1 
    250                pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    251             END DO 
    252          END DO 
    253          DO jk = jpk-2, 1, -1 
    254             DO jj = 2, jpjm1 
    255                DO ji = fs_2, fs_jpim1 
    256                   pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
    257                      &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    258                END DO 
    259             END DO 
    260          END DO 
     222         DO_2D( 0, 0, 0, 0 ) 
     223            pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     224         END_2D 
     225         DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 
     226            pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
     227               &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     228         END_3D 
    261229         !                                            ! ================= ! 
    262230      END DO                                          !  end tracer loop  ! 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/usrdef_hgr.F90

    r10074 r13710  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6163      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2] 
    6264      ! 
    63       INTEGER  ::   ji, jj   ! dummy loop indices 
     65      INTEGER  ::   ji, jj     ! dummy loop indices 
    6466      REAL(wp) ::   zphi0, zlam0, zbeta, zf0 
    65       REAL(wp) ::   zti, zui, ztj, zvj   ! local scalars 
     67      REAL(wp) ::   zti, ztj   ! local scalars 
    6668      !!------------------------------------------------------------------------------- 
    6769      ! 
    … …  
    7577      ! Position coordinates (in kilometers) 
    7678      !                          ========== 
    77       zlam0 = -REAL(NINT(jpiglo*rn_0xratio)-1, wp) * rn_dx 
    78       zphi0 = -REAL(NINT(jpjglo*rn_0yratio)-1, wp) * rn_dy  
     79      zlam0 = -REAL(NINT(Ni0glo*rn_0xratio)-1, wp) * rn_dx 
     80      zphi0 = -REAL(NINT(Nj0glo*rn_0yratio)-1, wp) * rn_dy  
    7981 
    8082#if defined key_agrif 
    … …  
    8890#endif 
    8991          
    90       DO jj = 1, jpj 
    91          DO ji = 1, jpi 
    92             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    93             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    94  
    95             plamt(ji,jj) = zlam0 + rn_dx * zti 
    96             plamu(ji,jj) = zlam0 + rn_dx * zui 
    97             plamv(ji,jj) = plamt(ji,jj)  
    98             plamf(ji,jj) = plamu(ji,jj)  
    99     
    100             pphit(ji,jj) = zphi0 + rn_dy * ztj 
    101             pphiv(ji,jj) = zphi0 + rn_dy * zvj 
    102             pphiu(ji,jj) = pphit(ji,jj)  
    103             pphif(ji,jj) = pphiv(ji,jj)  
    104          END DO 
    105       END DO 
     92      DO_2D( 1, 1, 1, 1 )          
     93         zti = REAL( mig0_oldcmp(ji) - 1, wp )   ! start at i=0 in the global grid without halos 
     94         ztj = REAL( mjg0_oldcmp(jj) - 1, wp )   ! start at j=0 in the global grid without halos 
     95          
     96         plamt(ji,jj) = zlam0 + rn_dx *   zti 
     97         plamu(ji,jj) = zlam0 + rn_dx * ( zti + 0.5_wp )  
     98         plamv(ji,jj) = plamt(ji,jj)  
     99         plamf(ji,jj) = plamu(ji,jj)  
     100          
     101         pphit(ji,jj) = zphi0 + rn_dy *   ztj 
     102         pphiv(ji,jj) = zphi0 + rn_dy * ( ztj + 0.5_wp )  
     103         pphiu(ji,jj) = pphit(ji,jj)  
     104         pphif(ji,jj) = pphiv(ji,jj)  
     105      END_2D 
    106106      !      
    107107      ! Horizontal scale factors (in meters) 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/usrdef_istate.F90

    r12489 r13710  
    2828   PUBLIC   usr_def_istate   ! called by istate.F90 
    2929 
     30   !! * Substitutions 
     31#  include "do_loop_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6466      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   ' 
    6567      ! 
    66       IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom) 
    6768      zjetx = ABS(rn_ujetszx)/2. 
    6869      zjety = ABS(rn_ujetszy)/2. 
    6970      ! 
     71      zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 
     72      ! 
    7073      SELECT CASE(nn_initcase) 
     74 
     75      CASE(-1)    ! stratif at rest 
     76 
     77         ! sea level: 
     78         pssh(:,:) = 0. 
     79         ! temperature: 
     80         pts(:,:,1,jp_tem) = 25. !!30._wp 
     81         pts(:,:,2:jpk,jp_tem) = 22. !!24._wp 
     82         ! salinity:   
     83         pts(:,:,:,jp_sal) = 35._wp 
     84         ! velocities: 
     85         pu(:,:,:) = 0. 
     86         pv(:,:,:) = 0. 
     87 
    7188      CASE(0)    ! rest 
    7289          
    … …  
    96113            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 
    97114            WHERE( ABS(gphit) <= zjety ) 
    98                pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
    99             ELSEWHERE 
    100                pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   & 
     115               pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
     116            ELSEWHERE 
     117               pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3   & 
    101118                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
    102119            END WHERE 
    … …  
    107124         pts(:,:,jpk,jp_sal) = 0. 
    108125         DO jk=1, jpkm1 
    109             pts(:,:,jk,jp_sal) = gphit(:,:) 
     126            WHERE( ABS(gphit) <= zjety ) 
     127!!$            WHERE( ABS(gphit) <= zjety*0.5 .AND. ABS(glamt) <= zjety*0.5 ) ! for a square of salt 
     128               pts(:,:,jk,jp_sal) = 35. 
     129            ELSEWHERE 
     130               pts(:,:,jk,jp_sal) = 30. 
     131            END WHERE                     
    110132         END DO 
    111133         ! velocities: 
    … …  
    132154            WHERE( ABS(gphit) <= zjety ) 
    133155               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   & 
    134                   &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
     156                  &        * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
    135157            ELSEWHERE 
    136158               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   & 
    137                   &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
     159                  &        * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
    138160            END WHERE 
    139161         END SELECT 
    … …  
    141163         pts(:,:,:,jp_tem) = 10._wp 
    142164         ! salinity:   
    143          pts(:,:,:,jp_sal) = 2. 
    144          DO jk=1, jpkm1 
    145             WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 2. + SIGN(1.,gphiv(:,:)) 
     165         pts(:,:,:,jp_sal) = 30. 
     166         DO jk=1, jpkm1 
     167            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 30. + SIGN(1.,gphiv(:,:)) 
    146168         END DO 
    147169         ! velocities: 
    … …  
    164186         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 
    165187         DO jl=1, jpnj 
    166             DO jj=nldj, nlej 
    167                DO ji=nldi, nlei 
    168                   pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
    169                END DO 
    170             END DO 
     188            DO_2D( 0, 0, 0, 0 ) 
     189               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
     190            END_2D 
    171191            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    172192         END DO 
    … …  
    176196         ! salinity:   
    177197         DO jk=1, jpkm1 
    178             pts(:,:,jk,jp_sal) = gphit(:,:) 
     198            pts(:,:,jk,jp_sal) = pssh(:,:) 
    179199         END DO 
    180200         ! velocities: 
    … …  
    183203      CASE(4)    ! geostrophic zonal pulse 
    184204    
    185          DO jj=1, jpj 
    186             DO ji=1, jpi 
    187                IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
    188                   zdu = rn_uzonal 
    189                ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
    190                   zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
    191                ELSE 
    192                   zdu = 0. 
    193                END IF 
    194                IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
    195                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
    196                   pu(ji,jj,:) = zdu 
    197                   pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
    198                ELSE 
    199                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
    200                   pu(ji,jj,:) = 0. 
    201                   pts(ji,jj,:,jp_sal) = 1. 
    202                END IF 
    203             END DO 
    204          END DO 
     205         DO_2D( 1, 1, 1, 1 ) 
     206            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
     207               zdu = rn_uzonal 
     208            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
     209               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
     210            ELSE 
     211               zdu = 0. 
     212            END IF 
     213            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
     214               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
     215               pu(ji,jj,:) = zdu 
     216               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
     217            ELSE 
     218               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
     219               pu(ji,jj,:) = 0. 
     220               pts(ji,jj,:,jp_sal) = 1. 
     221            END IF 
     222         END_2D 
    205223          
    206224         ! temperature: 
    207225         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)         
    208226         pv(:,:,:) = 0. 
    209           
    210227          
    211228       CASE(5)    ! vortex 
    … …  
    213230         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 
    214231         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic 
    215          zlambda = SQRT(2._wp)*rn_lambda       ! Horizontal scale in meters  
     232         zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters  
    216233         zn2 = 3.e-3**2 
    217234         zH = 0.5_wp * 5000._wp 
    … …  
    220237         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 
    221238         ! 
    222          DO jj=1, jpj 
    223             DO ji=1, jpi 
    224                zx = glamt(ji,jj) * 1.e3 
    225                zy = gphit(ji,jj) * 1.e3 
    226                ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
    227                zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
    228                ! Sea level: 
    229                pssh(ji,jj) = 0. 
    230                DO jl=1,5 
    231                   zdt = pssh(ji,jj) 
    232                   zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
    233                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    234                   pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
    235                END DO 
    236                ! temperature: 
    237                DO jk=1,jpk 
    238                   zdt =  pdept(ji,jj,jk)  
    239                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
    240                   IF (zdt < zH) THEN 
    241                      zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
    242                      zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    243                   ENDIF 
    244                   !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    245                   pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    246                END DO 
     239         DO_2D( 1, 1, 1, 1 ) 
     240            zx = glamt(ji,jj) * 1.e3 
     241            zy = gphit(ji,jj) * 1.e3 
     242            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
     243            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
     244            ! Sea level: 
     245            pssh(ji,jj) = 0. 
     246            DO jl=1,5 
     247               zdt = pssh(ji,jj) 
     248               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
     249               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     250               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
    247251            END DO 
    248          END DO 
     252            ! temperature: 
     253            DO jk=1,jpk 
     254               zdt =  pdept(ji,jj,jk)  
     255               zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
     256               IF (zdt < zH) THEN 
     257                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
     258                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     259               ENDIF 
     260               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     261               pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     262            END DO 
     263         END_2D 
    249264         ! 
    250265         ! salinity:   
    … …  
    253268         ! velocities: 
    254269         za = 2._wp * zP0 / zlambda**2 
    255          DO jj=1, jpj 
    256             DO ji=1, jpim1 
    257                zx = glamu(ji,jj) * 1.e3 
    258                zy = gphiu(ji,jj) * 1.e3 
    259                DO jk=1, jpk 
    260                   zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
    261                   IF (zdu < zH) THEN 
    262                      zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
    263                      zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
    264                      pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
    265                   ELSE 
    266                      pu(ji,jj,jk) = 0._wp 
    267                   ENDIF 
    268                END DO 
     270         DO_2D( 0, 0, 0, 0 ) 
     271            zx = glamu(ji,jj) * 1.e3 
     272            zy = gphiu(ji,jj) * 1.e3 
     273            DO jk=1, jpk 
     274               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
     275               IF (zdu < zH) THEN 
     276                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
     277                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
     278                  pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
     279               ELSE 
     280                  pu(ji,jj,jk) = 0._wp 
     281               ENDIF 
    269282            END DO 
    270          END DO 
    271          ! 
    272          DO jj=1, jpjm1 
    273             DO ji=1, jpi 
    274                zx = glamv(ji,jj) * 1.e3 
    275                zy = gphiv(ji,jj) * 1.e3 
    276                DO jk=1, jpk 
    277                   zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
    278                   IF (zdv < zH) THEN 
    279                      zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
    280                      zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
    281                      pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
    282                   ELSE 
    283                      pv(ji,jj,jk) = 0._wp 
    284                   ENDIF 
    285                END DO 
     283         END_2D 
     284         ! 
     285         DO_2D( 0, 0, 0, 0 ) 
     286            zx = glamv(ji,jj) * 1.e3 
     287            zy = gphiv(ji,jj) * 1.e3 
     288            DO jk=1, jpk 
     289               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
     290               IF (zdv < zH) THEN 
     291                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
     292                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
     293                  pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
     294               ELSE 
     295                  pv(ji,jj,jk) = 0._wp 
     296               ENDIF 
    286297            END DO 
    287          END DO 
     298         END_2D 
    288299         !             
    289300      END SELECT 
    290  
     301       
    291302      IF (ln_sshnoise) THEN 
     303         CALL RANDOM_SEED() 
    292304         CALL RANDOM_NUMBER(zrandom) 
    293305         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 ) 
    294306      END IF 
    295307      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    296       CALL lbc_lnk(  'usrdef_istate', pts, 'T',  1. ) 
    297       CALL lbc_lnk(   'usrdef_istate', pu, 'U', -1. ) 
    298       CALL lbc_lnk(   'usrdef_istate', pv, 'V', -1. ) 
     308      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. ) 
     309      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 
    299310 
    300311   END SUBROUTINE usr_def_istate 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/usrdef_nam.F90

    r12377 r13710  
    1414   !!   usr_def_hgr   : initialize the horizontal mesh  
    1515   !!---------------------------------------------------------------------- 
    16    USE dom_oce  , ONLY: nimpp , njmpp            ! i- & j-indices of the local domain 
     16   USE dom_oce 
    1717   USE par_oce        ! ocean space and time domain 
    1818   USE phycst         ! physical constants 
    … …  
    5050   LOGICAL , PUBLIC ::   ln_sshnoise=.false. ! add random noise on initial ssh 
    5151   REAL(wp), PUBLIC ::   rn_lambda  = 50.    ! gaussian lambda 
     52   INTEGER , PUBLIC ::   nn_perio   =    0   ! periodicity of the channel (0=closed, 1=E-W) 
    5253 
    5354   !!---------------------------------------------------------------------- 
    … …  
    7980      !! 
    8081      NAMELIST/namusr_def/  rn_domszx, rn_domszy, rn_domszz, rn_dx, rn_dy, rn_dz, rn_0xratio, rn_0yratio   & 
    81          &                 , nn_fcase, rn_ppgphi0, rn_vtxmax, rn_uzonal, rn_ujetszx, rn_ujetszy   & 
    82          &                 , rn_u10, rn_windszx, rn_windszy, rn_uofac   & 
    83          &                 , nn_botcase, nn_initcase, ln_sshnoise, rn_lambda 
     82         &                 , nn_fcase, rn_ppgphi0, rn_u10, rn_windszx, rn_windszy & !!, rn_uofac   & 
     83         &                 , rn_vtxmax, rn_uzonal, rn_ujetszx, rn_ujetszy  & 
     84         &                 , nn_botcase, nn_initcase, ln_sshnoise, rn_lambda, nn_perio 
    8485      !!---------------------------------------------------------------------- 
    8586      ! 
    … …  
    106107      kk_cfg = INT( rn_dx ) 
    107108      ! 
    108       ! Global Domain size:  EW_CANAL global domain is  1800 km x 1800 Km x 5000 m 
    109       kpi = NINT( rn_domszx / rn_dx ) + 1 
    110       kpj = NINT( rn_domszy / rn_dy ) + 3 
    111       kpk = NINT( rn_domszz / rn_dz ) + 1 
    112 #if defined key_agrif 
    113       IF( .NOT. Agrif_Root() ) THEN 
    114          kpi  = nbcellsx + 2 + 2*nbghostcells 
    115          kpj  = nbcellsy + 2 + 2*nbghostcells 
     109      IF( Agrif_Root() ) THEN        ! Global Domain size:  EW_CANAL global domain is  1800 km x 1800 Km x 5000 m 
     110         kpi = NINT( rn_domszx / rn_dx ) + 1 
     111         kpj = NINT( rn_domszy / rn_dy ) + 3 
     112      ELSE                           ! Global Domain size: add nbghostcells + 1 "land" point on each side 
     113         kpi  = nbcellsx + nbghostcells_x   + nbghostcells_x   + 2 
     114         kpj  = nbcellsy + nbghostcells_y_s + nbghostcells_y_n + 2 
    116115      ENDIF 
    117 #endif 
     116      kpk = MAX( 2, NINT( rn_domszz / rn_dz ) + 1 ) 
    118117      ! 
    119118      zh  = (kpk-1)*rn_dz 
    … …  
    150149         WRITE(numout,*) '      add random noise on initial ssh   ln_sshnoise= ', ln_sshnoise 
    151150         WRITE(numout,*) '      Gaussian lambda parameter          rn_lambda = ', rn_lambda 
    152          WRITE(numout,*) '   ' 
    153          WRITE(numout,*) '   Lateral boundary condition of the global domain' 
    154          WRITE(numout,*) '      EW_CANAL : closed basin               jperio = ', kperio 
     151         WRITE(numout,*) '      Periodicity of the basin            nn_perio = ', nn_perio 
    155152      ENDIF 
     153      !                             ! Set the lateral boundary condition of the global domain 
     154      kperio = nn_perio                    ! EW_CANAL configuration : closed basin 
    156155      ! 
    157156   END SUBROUTINE usr_def_nam 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/usrdef_sbc.F90

    r12377 r13710  
    1717   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1818   USE phycst          ! physical constants 
    19    USE usrdef_nam, ONLY : rn_u10, rn_uofac, rn_windszy  
     19   USE usrdef_nam, ONLY : rn_u10, rn_uofac, rn_windszy, rn_windszx  
    2020   ! 
    2121   USE in_out_manager  ! I/O manager 
    … …  
    3838CONTAINS 
    3939 
    40    SUBROUTINE usrdef_sbc_oce( kt, Kmm, Kbb ) 
     40   SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 
    4141      !!--------------------------------------------------------------------- 
    4242      !!                    ***  ROUTINE usr_def_sbc  *** 
    … …  
    5353      !!---------------------------------------------------------------------- 
    5454      INTEGER, INTENT(in) ::   kt        ! ocean time step 
    55       INTEGER, INTENT(in) ::   Kbb, Kmm  ! ocean time index 
     55      INTEGER, INTENT(in) ::   Kbb       ! ocean time index 
    5656      INTEGER  ::   ji, jj               ! dummy loop indices 
    5757      REAL(wp) :: zrhoair = 1.22     ! approximate air density [Kg/m3] 
    … …  
    6969         ! 
    7070         utau(:,:) = 0._wp 
    71          IF( rn_u10 /= 0. .AND. rn_windszy > 0. ) THEN 
    72             WHERE( ABS(gphit) <= rn_windszy/2. ) utau(:,:) = zrhocd * rn_u10 * rn_u10 
    73          ENDIF 
    7471         vtau(:,:) = 0._wp 
    7572         taum(:,:) = 0._wp 
    … …  
    8380      ENDIF 
    8481 
     82      IF( rn_u10 /= 0. .AND. rn_windszy > 0. ) THEN 
     83         IF( nyear == 1 .AND. nmonth == 1 .AND. nday <= 10 ) THEN 
     84            WHERE( ABS(gphit) <= rn_windszy/2. .AND. ABS(glamt) <= rn_windszx/2. ) utau(:,:) = zrhocd * rn_u10 * rn_u10 
     85         ELSE 
     86            utau(:,:) = 0. 
     87         ENDIF 
     88      ENDIF 
     89 
    8590      IF( rn_uofac /= 0. ) THEN 
    8691          
    8792         WHERE( ABS(gphit) <= rn_windszy/2. ) 
    88             zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kmm) 
     93            zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kbb) 
    8994         ELSEWHERE 
    90             zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kmm) 
     95            zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kbb) 
    9196         END WHERE 
    9297         utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    9398 
    94          zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kmm) 
     99         zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kbb) 
    95100         vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    96101 

  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/tests/CANAL/MY_SRC/usrdef_zgr.F90

    r12377 r13710  
    197197         zmaxlam = MAXVAL(glamt) 
    198198         CALL mpp_max( 'usrdef_zgr', zmaxlam )                 ! max over the global domain 
    199          zscl = rpi / zmaxlam 
    200          z2d(:,:) = 0.5 * ( 1. - COS( glamt(:,:) * zscl ) ) 
    201          z2d(:,:) = REAL(jpkm1 - NINT( 0.75 * REAL(jpkm1,wp) * z2d(:,:) ), wp) 
     199         zscl = 0.5 * rpi / zmaxlam 
     200         z2d(:,:) = COS( glamt(:,:) * zscl ) 
     201         z2d(:,:) = REAL(jpkm1 - NINT( 0.5 * REAL(jpkm1,wp) * z2d(:,:) ), wp) 
    202202      END SELECT 
    203203      ! 
    204204      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    205205      ! 
    206       k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
     206      k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere 
    207207      ! 
    208208      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere 
Note: See TracChangeset for help on using the changeset viewer.