Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/context_nemo.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/context_nemo.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/context_nemo.xml (revision 10118)
@@ -0,0 +1,35 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/domain_def_nemo.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/domain_def_nemo.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/domain_def_nemo.xml (revision 10118)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/domain_def_nemo.xml
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/field_def_nemo-oce.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/field_def_nemo-oce.xml (revision 10118)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-oce.xml
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/file_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/file_def_nemo-oce.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/file_def_nemo-oce.xml (revision 10118)
@@ -0,0 +1,36 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/grid_def_nemo.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/grid_def_nemo.xml (revision 10118)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/grid_def_nemo.xml
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/iodef.xml
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/iodef.xml (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/iodef.xml (revision 10118)
@@ -0,0 +1,25 @@
+
+
+
+
+
+
+
+
+
+
+
+ 10
+ true
+ false
+ oceanx
+
+
+
+
+
+
+
+
+
+
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_cfg (revision 10118)
@@ -0,0 +1,350 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OPA Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! COMODO-IW configuration !!
+!!======================================================================
+!! *** Domain & Run management namelists *** !!
+!! !!
+!! namrun parameters of the run
+!! namdom space and time domain
+!! namcfg parameters of the configuration (default: user defined GYRE)
+!! namwad Wetting and drying (default: OFF)
+!! namtsd data: temperature & salinity (default: OFF)
+!! namcrs coarsened grid (for outputs and/or TOP) (ln_crs =T)
+!! namc1d 1D configuration options ("key_c1d")
+!! namc1d_dyndmp 1D newtonian damping applied on currents ("key_c1d")
+!! namc1d_uvd 1D data (currents) ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namcfg ! parameters of the configuration (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+ ln_read_cfg = .true. ! (=T) read the domain configuration file
+/
+!-----------------------------------------------------------------------
+&namrun ! parameters of the run
+!-----------------------------------------------------------------------
+ nn_no = 0 ! job number (no more used...)
+ cn_exp = "CO-IW" ! experience name
+ nn_it000 = 1 ! first time step
+ nn_itend = 21600 ! last time step
+ nn_istate = 0 ! output the initial state (1) or not (0)
+ nn_stock = 99999 ! frequency of creation of a restart file (modulo referenced to 1)
+ nn_write = 99999 ! frequency of write in the output file (modulo referenced to nn_it000)
+/
+!-----------------------------------------------------------------------
+&namdom ! time and space domain
+!-----------------------------------------------------------------------
+ ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time
+ rn_rdt = 120. ! time step for the dynamics (and tracer if nn_acc=0)
+ rn_atfp = 0.05 ! asselin time filter parameter
+/
+!-----------------------------------------------------------------------
+&namcfg ! parameters of the configuration (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+/
+!!======================================================================
+!! *** Surface Boundary Condition namelists *** !!
+!! !!
+!! namsbc surface boundary condition manager (default: NO selection)
+!! namsbc_flx flux formulation (ln_flx =T)
+!! namsbc_blk Bulk formulae formulation (ln_blk =T)
+!! namsbc_cpl CouPLed formulation ("key_oasis3" )
+!! namsbc_sas Stand-Alone Surface module (SAS_SRC only)
+!! namsbc_iif Ice-IF: use observed ice cover (nn_ice = 1 )
+!! namtra_qsr penetrative solar radiation (ln_traqsr =T)
+!! namsbc_ssr sea surface restoring term (for T and/or S) (ln_ssr =T)
+!! namsbc_rnf river runoffs (ln_rnf =T)
+!! namsbc_apr Atmospheric Pressure (ln_apr_dyn =T)
+!! namsbc_isf ice shelf melting/freezing (ln_isfcav =T : read (ln_read_cfg=T) or set or usr_def_zgr )
+!! namsbc_iscpl coupling option between land ice model and ocean (ln_isfcav =T)
+!! namsbc_wave external fields from wave model (ln_wave =T)
+!! namberg iceberg floats (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc ! Surface Boundary Condition manager (default: NO selection)
+!-----------------------------------------------------------------------
+ nn_fsbc = 1 ! frequency of surface boundary condition computation
+ ! (also = the frequency of sea-ice & iceberg model call)
+ ln_usr = .true. ! user defined formulation (T => check usrdef_sbc)
+ ln_blk = .false. ! Bulk formulation (T => fill namsbc_blk )
+ nn_ice = 0 ! =0 no ice boundary condition
+ ln_traqsr = .false. ! Light penetration in the ocean (T => fill namtra_qsr )
+ ln_rnf = .false. ! runoffs (T => fill namsbc_rnf)
+ ln_ssr = .false. ! Sea Surface Restoring on T and/or S (T => fill namsbc_ssr)
+ nn_fwb = 0 ! FreshWater Budget: =0 unchecked
+/
+!!======================================================================
+!! *** Lateral boundary condition *** !!
+!! !!
+!! namlbc lateral momentum boundary condition (default: NO selection)
+!! namagrif agrif nested grid ( read by child model only ) ("key_agrif")
+!! nam_tide Tidal forcing (default: OFF)
+!! nambdy Unstructured open boundaries (default: OFF)
+!! nambdy_dta Unstructured open boundaries - external data (see nambdy)
+!! nambdy_tide tidal forcing at open boundaries (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc ! lateral momentum boundary condition (default: NO selection)
+!-----------------------------------------------------------------------
+ rn_shlat = 0. ! free slip
+/
+!!======================================================================
+!! *** Top/Bottom boundary condition *** !!
+!! !!
+!! namdrg top/bottom drag coefficient (default: NO selection)
+!! namdrg_top top friction (ln_OFF=F & ln_isfcav=T)
+!! namdrg_bot bottom friction (ln_OFF=F)
+!! nambbc bottom temperature boundary condition (default: OFF)
+!! nambbl bottom boundary layer scheme (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdrg ! top/bottom drag coefficient (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_OFF = .true. ! free-slip : Cd = 0
+/
+!-----------------------------------------------------------------------
+&nam_tide ! tide parameters (default: OFF)
+!-----------------------------------------------------------------------
+ ln_tide = .true. ! Activate tides
+ ln_tide_pot = .false. ! use tidal potential forcing
+ ln_tide_ramp = .true. ! Use linear ramp for tides at startup
+ rdttideramp = 2. ! ramp duration in days
+ clname(1) = 'S2' ! name of constituent - all tidal components must be set in namelist_cfg
+/
+!-----------------------------------------------------------------------
+&nambdy ! unstructured open boundaries (default: OFF)
+!-----------------------------------------------------------------------
+ ln_bdy = .true. ! Use unstructured open boundaries
+ nb_bdy = 2 ! number of open boundary sets
+ ln_coords_file = .false. ! =T : read bdy coordinates from file
+ cn_dyn2d = 'flather','flather' !
+ nn_dyn2d_dta = 2,2 ! = 0, bdy data are equal to the initial state
+ ! ! = 1, bdy data are read in 'bdydata .nc' files
+ ! ! = 2, use tidal harmonic forcing data from files
+ ! ! = 3, use external data AND tidal harmonic forcing
+ cn_dyn3d = 'specified','specified' !
+ nn_dyn3d_dta = 0,0 ! = 0, bdy data are equal to the initial state
+ ! ! = 1, bdy data are read in 'bdydata .nc' files
+ cn_tra = 'specified','specified' !
+ nn_tra_dta = 0,0 ! = 0, bdy data are equal to the initial state
+ ! ! = 1, bdy data are read in 'bdydata .nc' files
+ !
+ ln_tra_dmp =.true.,.true. ! open boudaries conditions for tracers
+ ln_dyn3d_dmp =.true.,.true. ! open boundary condition for baroclinic velocities
+ rn_time_dmp = 0.05,0.05 ! Damping time scale in days
+ nn_rimwidth = 155,155 ! width of the relaxation zone
+/
+!-----------------------------------------------------------------------
+&nambdy_index ! structured open boundaries definition ("key_bdy")
+!-----------------------------------------------------------------------
+ ctypebdy ='W' ! Open boundary type (W,E,S or N)
+ nbdyind = -1 ! indice of velocity row or column
+ ! if ==-1, set obc at the domain boundary
+ ! , discard start and end indices
+ nbdybeg = -99 ! indice of segment start
+ nbdyend = -99 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index ! structured open boundaries definition ("key_bdy")
+!-----------------------------------------------------------------------
+ ctypebdy ='E' ! Open boundary type (W,E,S or N)
+ nbdyind = -1 ! indice of velocity row or column
+ ! if ==-1, set obc at the domain boundary
+ ! , discard start and end indices
+ nbdybeg = -99 ! indice of segment start
+ nbdyend = -99 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_tide ! tidal forcing at open boundaries (default: OFF)
+!-----------------------------------------------------------------------
+ filtide = 'tide' ! file name root of tidal forcing files
+ ln_bdytide_2ddta = .true. !
+ ln_bdytide_conj = .true. !
+/
+!-----------------------------------------------------------------------
+&nambdy_tide ! tidal forcing at open boundaries (default: OFF)
+!-----------------------------------------------------------------------
+ filtide = 'tide' ! file name root of tidal forcing files
+ ln_bdytide_2ddta = .true. !
+ ln_bdytide_conj = .true. !
+/
+!!======================================================================
+!! Tracer (T-S) namelists !!
+!! !!
+!! nameos equation of state (default: NO selection)
+!! namtra_adv advection scheme (default: NO selection)
+!! namtra_ldf lateral diffusion scheme (default: NO selection)
+!! namtra_mle mixed layer eddy param. (Fox-Kemper param.) (default: OFF)
+!! namtra_eiv eddy induced velocity param. (default: OFF)
+!! namtra_dmp T & S newtonian damping (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos ! ocean Equation Of Seawater (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_seos = .true. ! = Use simplified equation of state (S-EOS)
+ ! ! rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+ rn_a0 = 0.16550 ! thermal expension coefficient (for simplified equation of state)
+ rn_b0 = 0. ! saline expension coefficient (for simplified equation of state)
+ rn_lambda1 = 0. ! cabbeling coeff in T^2 (=0 for linear eos)
+ rn_lambda2 = 0. ! cabbeling coeff in S^2 (=0 for linear eos)
+ rn_mu1 = 0. ! thermobaric coeff. in T (=0 for linear eos)
+ rn_mu2 = 0. ! thermobaric coeff. in S (=0 for linear eos)
+ rn_nu = 0. ! cabbeling coeff in T*S (=0 for linear eos)
+/
+!-----------------------------------------------------------------------
+&namtra_adv ! advection scheme for tracer (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_traadv_fct = .true. ! FCT scheme
+ nn_fct_h = 2 ! =2/4, horizontal 2nd / 4th order
+ nn_fct_v = 2 ! =2/4, vertical 2nd / COMPACT 4th order
+ ln_traadv_qck = .false. ! QUICKEST scheme
+/
+!-----------------------------------------------------------------------
+&namtra_ldf ! lateral diffusion scheme for tracers (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_traldf_OFF = .true. ! No explicit diffusion
+/
+!!======================================================================
+!! *** Dynamics namelists *** !!
+!! !!
+!! nam_vvl vertical coordinate options (default: z-star)
+!! namdyn_adv formulation of the momentum advection (default: NO selection)
+!! namdyn_vor advection scheme (default: NO selection)
+!! namdyn_hpg hydrostatic pressure gradient (default: NO selection)
+!! namdyn_spg surface pressure gradient (default: NO selection)
+!! namdyn_ldf lateral diffusion scheme (default: NO selection)
+!! namdta_dyn offline TOP: dynamics read in files (OFF_SRC only)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nam_vvl ! vertical coordinate options (default: z-star)
+!-----------------------------------------------------------------------
+ ln_vvl_zstar = .false. ! zstar vertical coordinate
+ ln_vvl_ztilde = .true. ! ztilde vertical coordinate: only high frequency variations
+ ln_vvl_layer = .false. ! full layer vertical coordinate
+ ln_vvl_ztilde_as_zstar = .false. ! ztilde vertical coordinate emulating zstar
+ ln_vvl_zstar_at_eqtor = .false. ! revert to zstart at equator
+ ln_vvl_zstar_on_shelf = .false. ! revert to zstar on shelves
+ ln_vvl_adv_cn2 = .false. ! Centred thickness advection scheme
+ ln_vvl_adv_fct = .true. ! FCT thickness advection scheme
+ ln_vvl_lap = .false. ! Laplacian thickness diffusion
+ rn_ahe3_lap = 30.e0
+ ln_vvl_blp = .false. ! Bilaplacian thickness diffusion
+ rn_ahe3_blp = -2.e8
+ rn_rst_e3t = 30.e0 ! ztilde to zstar restoration timescale [days]
+ rn_lf_cutoff = 5.0e0 ! cutoff frequency for low-pass filter [days]
+ ln_vvl_regrid = .true. ! Vertical regriding to ensure layer separation
+ ln_vvl_dbg = .false. ! debug prints (T/F)
+ ln_vvl_ramp = .false. ! linear ramp at startup
+ rn_day_ramp = 3.e0 ! ramp duration [days]
+/
+!-----------------------------------------------------------------------
+&namdyn_adv ! formulation of the momentum advection (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_dynadv_ubs = .true. ! flux form - 3rd order UBS scheme
+/
+!-----------------------------------------------------------------------
+&namdyn_vor ! Vorticity / Coriolis scheme (default: OFF)
+!-----------------------------------------------------------------------
+ ln_dynvor_een = .true. ! energy & enstrophy scheme
+ nn_een_e3f = 1 ! e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg ! Hydrostatic pressure gradient option (default: NO selection)
+!-----------------------------------------------------------------------
+ ln_hpg_sco = .true. ! s-coordinate (standard jacobian formulation)
+/
+!-----------------------------------------------------------------------
+&namdyn_spg ! surface pressure gradient (default: OFF)
+!-----------------------------------------------------------------------
+ ln_dynspg_exp = .false. ! explicit free surface
+ ln_dynspg_ts = .true. ! split-explicit free surface
+ ln_bt_fw = .false. ! Forward integration of barotropic Eqs.
+ ln_bt_av = .true. ! Time filtering of barotropic variables
+ nn_bt_flt = 1 ! Time filter choice = 0 None
+ ! ! = 1 Boxcar over nn_baro sub-steps
+ ! ! = 2 Boxcar over 2*nn_baro " "
+ ln_bt_auto = .false. ! Number of sub-step defined from:
+ nn_baro = 48 ! =F : the number of sub-step in rn_rdt seconds
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf ! lateral diffusion on momentum (default: NO selection)
+!-----------------------------------------------------------------------
+ ! ! Type of the operator :
+ ln_dynldf_OFF = .true. ! No operator (i.e. no explicit diffusion)
+/
+!!======================================================================
+!! vertical physics namelists !!
+!! !!
+!! namzdf vertical physics manager (default: NO selection)
+!! namzdf_ric richardson number vertical mixing (ln_zdfric=T)
+!! namzdf_tke TKE vertical mixing (ln_zdftke=T)
+!! namzdf_gls GLS vertical mixing (ln_zdfgls=T)
+!! namzdf_osm OSM vertical diffusion (ln_zdfosm=T)
+!! namzdf_iwm tidal mixing parameterization (ln_zdfiwm=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namzdf ! vertical physics manager (default: NO selection)
+!-----------------------------------------------------------------------
+ ! ! type of vertical closure
+ ln_zdfcst = .true. ! constant mixing
+ !
+ ! ! convection
+ ln_zdfevd = .false. ! enhanced vertical diffusion
+ ln_zdfnpc = .false. ! Non-Penetrative Convective algorithm
+ !
+ ! ! coefficients
+ rn_avm0 = 0.e0 ! vertical eddy viscosity [m2/s] (background Kz if ln_zdfcst=F)
+ rn_avt0 = 0.e0 ! vertical eddy diffusivity [m2/s] (background Kz if ln_zdfcst=F)
+ nn_avb = 0 ! profile for background avt & avm (=1) or not (=0)
+ nn_havtb = 0 ! horizontal shape for avtb (=1) or not (=0)
+/
+!!======================================================================
+!! *** Diagnostics namelists *** !!
+!! !!
+!! namtrd dynamics and/or tracer trends (default: OFF)
+!! namptr Poleward Transport Diagnostics (default: OFF)
+!! namhsb Heat and salt budgets (default: OFF)
+!! namdiu Cool skin and warm layer models (default: OFF)
+!! namdiu Cool skin and warm layer models (default: OFF)
+!! namflo float parameters ("key_float")
+!! nam_diaharm Harmonic analysis of tidal constituents ("key_diaharm")
+!! namdct transports through some sections ("key_diadct")
+!! nam_diatmb Top Middle Bottom Output (default: OFF)
+!! nam_dia25h 25h Mean Output (default: OFF)
+!! namnc4 netcdf4 chunking and compression settings ("key_netcdf4")
+!!======================================================================
+!
+!!======================================================================
+!! *** Observation & Assimilation *** !!
+!! !!
+!! namobs observation and model comparison (default: OFF)
+!! nam_asminc assimilation increments ('key_asminc')
+!!======================================================================
+!
+!!======================================================================
+!! *** Miscellaneous namelists *** !!
+!! !!
+!! nammpp Massively Parallel Processing ("key_mpp_mpi")
+!! namctl Control prints (default: OFF)
+!! namsto Stochastic parametrization of EOS (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp ! Massively Parallel Processing ("key_mpp_mpi")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl ! Control prints (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsto ! Stochastic parametrization of EOS (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_ref
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_ref (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/EXPREF/namelist_ref (revision 10118)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ref
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/bdyini.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/bdyini.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/bdyini.F90 (revision 10118)
@@ -0,0 +1,1779 @@
+MODULE bdyini
+ !!======================================================================
+ !! *** MODULE bdyini ***
+ !! Unstructured open boundaries : initialisation
+ !!======================================================================
+ !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
+ !! - ! 2007-01 (D. Storkey) Update to use IOM module
+ !! - ! 2007-01 (D. Storkey) Tidal forcing
+ !! 3.0 ! 2008-04 (NEMO team) add in the reference version
+ !! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations
+ !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
+ !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge
+ !! 3.4 ! 2012 (J. Chanut) straight open boundary case update
+ !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) optimization of BDY communications
+ !! 3.7 ! 2016 (T. Lovato) Remove bdy macro, call here init for dta and tides
+ !!----------------------------------------------------------------------
+ !! bdy_init : Initialization of unstructured open boundaries
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers variables
+ USE dom_oce ! ocean space and time domain
+ USE bdy_oce ! unstructured open boundary conditions
+ USE bdydta ! open boundary cond. setting (bdy_dta_init routine)
+ USE bdytides ! open boundary cond. setting (bdytide_init routine)
+ USE sbctide ! Tidal forcing or not
+ USE phycst , ONLY: rday
+ !
+ USE in_out_manager ! I/O units
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_mpp ! for mpp_sum
+ USE iom ! I/O
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC bdy_init ! routine called in nemo_init
+
+ INTEGER, PARAMETER :: jp_nseg = 100 !
+ INTEGER, PARAMETER :: nrimmax = 155 ! maximum rimwidth in structured
+ ! open boundary data files
+ ! Straight open boundary segment parameters:
+ INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs
+ INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge !
+ INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw !
+ INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn !
+ INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs !
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id: bdyini.F90 9807 2018-06-15 22:51:16Z lovato $
+ !! Software governed by the CeCILL licence (./LICENSE)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE bdy_init
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE bdy_init ***
+ !!
+ !! ** Purpose : Initialization of the dynamics and tracer fields with
+ !! unstructured open boundaries.
+ !!
+ !! ** Method : Read initialization arrays (mask, indices) to identify
+ !! an unstructured open boundary
+ !!
+ !! ** Input : bdy_init.nc, input file for unstructured open boundaries
+ !!----------------------------------------------------------------------
+ NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
+ & ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
+ & cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
+ & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
+ & cn_ice, nn_ice_dta, &
+ & rn_ice_tem, rn_ice_sal, rn_ice_age, &
+ & ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+ !
+ INTEGER :: ios ! Local integer output status for namelist read
+ !!----------------------------------------------------------------------
+
+ ! ------------------------
+ ! Read namelist parameters
+ ! ------------------------
+ REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries
+ READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
+ REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries
+ READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
+902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, nambdy )
+
+ IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE. ! forced for Agrif children
+
+ ! -----------------------------------------
+ ! unstructured open boundaries use control
+ ! -----------------------------------------
+ IF ( ln_bdy ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~'
+ !
+ ! Open boundaries definition (arrays and masks)
+ CALL bdy_segs
+ !
+ ! Open boundaries initialisation of external data arrays
+ CALL bdy_dta_init
+ !
+ ! Open boundaries initialisation of tidal harmonic forcing
+ IF( ln_tide ) CALL bdytide_init
+ !
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~'
+ !
+ ENDIF
+ !
+ END SUBROUTINE bdy_init
+
+
+ SUBROUTINE bdy_segs
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE bdy_init ***
+ !!
+ !! ** Purpose : Definition of unstructured open boundaries.
+ !!
+ !! ** Method : Read initialization arrays (mask, indices) to identify
+ !! an unstructured open boundary
+ !!
+ !! ** Input : bdy_init.nc, input file for unstructured open boundaries
+ !!----------------------------------------------------------------------
+ INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
+ INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers
+ INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - -
+ INTEGER :: igrd_start, igrd_end, jpbdta ! - -
+ INTEGER :: jpbdtau, jpbdtas ! - -
+ INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - -
+ INTEGER :: i_offset, j_offset ! - -
+ INTEGER , POINTER :: nbi, nbj, nbr ! short cuts
+ REAL(wp), POINTER :: flagu, flagv ! - -
+ REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields
+ REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
+ INTEGER, DIMENSION (2) :: kdimsz
+ INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays
+ INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
+ INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
+ CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
+ INTEGER :: com_east, com_west, com_south, com_north ! Flags for boundaries sending
+ INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving
+ INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates
+ REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat)
+ !!
+ CHARACTER(LEN=1) :: ctypebdy ! - -
+ INTEGER :: nbdyind, nbdybeg, nbdyend
+ !!
+ NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
+ INTEGER :: ios ! Local integer output status for namelist read
+ !!----------------------------------------------------------------------
+ !
+ cgrid = (/'t','u','v'/)
+
+ ! -----------------------------------------
+ ! Check and write out namelist parameters
+ ! -----------------------------------------
+! IF( jperio /= 0 ) CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,', &
+! & ' and general open boundary condition are not compatible' )
+
+ IF( nb_bdy == 0 ) THEN
+ IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
+ ELSE
+ IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
+ ENDIF
+
+ DO ib_bdy = 1,nb_bdy
+ IF(lwp) WRITE(numout,*) ' '
+ IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'
+
+ IF( ln_coords_file(ib_bdy) ) THEN
+ IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
+ ELSE
+ IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
+ ENDIF
+ IF(lwp) WRITE(numout,*)
+
+ IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '
+ SELECT CASE( cn_dyn2d(ib_bdy) )
+ CASE( 'none' )
+ IF(lwp) WRITE(numout,*) ' no open boundary condition'
+ dta_bdy(ib_bdy)%ll_ssh = .false.
+ dta_bdy(ib_bdy)%ll_u2d = .false.
+ dta_bdy(ib_bdy)%ll_v2d = .false.
+ CASE( 'frs' )
+ IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
+ dta_bdy(ib_bdy)%ll_ssh = .false.
+ dta_bdy(ib_bdy)%ll_u2d = .true.
+ dta_bdy(ib_bdy)%ll_v2d = .true.
+ CASE( 'flather' )
+ IF(lwp) WRITE(numout,*) ' Flather radiation condition'
+ dta_bdy(ib_bdy)%ll_ssh = .true.
+ dta_bdy(ib_bdy)%ll_u2d = .true.
+ dta_bdy(ib_bdy)%ll_v2d = .true.
+ CASE( 'orlanski' )
+ IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_ssh = .false.
+ dta_bdy(ib_bdy)%ll_u2d = .true.
+ dta_bdy(ib_bdy)%ll_v2d = .true.
+ CASE( 'orlanski_npo' )
+ IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_ssh = .false.
+ dta_bdy(ib_bdy)%ll_u2d = .true.
+ dta_bdy(ib_bdy)%ll_v2d = .true.
+ CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
+ END SELECT
+ IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
+ SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !
+ CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
+ CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
+ CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file'
+ CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'
+ CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
+ END SELECT
+ IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN
+ CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
+ ENDIF
+ ENDIF
+ IF(lwp) WRITE(numout,*)
+
+ IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
+ SELECT CASE( cn_dyn3d(ib_bdy) )
+ CASE('none')
+ IF(lwp) WRITE(numout,*) ' no open boundary condition'
+ dta_bdy(ib_bdy)%ll_u3d = .false.
+ dta_bdy(ib_bdy)%ll_v3d = .false.
+ CASE('frs')
+ IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
+ dta_bdy(ib_bdy)%ll_u3d = .true.
+ dta_bdy(ib_bdy)%ll_v3d = .true.
+ CASE('specified')
+ IF(lwp) WRITE(numout,*) ' Specified value'
+ dta_bdy(ib_bdy)%ll_u3d = .true.
+ dta_bdy(ib_bdy)%ll_v3d = .true.
+ CASE('neumann')
+ IF(lwp) WRITE(numout,*) ' Neumann conditions'
+ dta_bdy(ib_bdy)%ll_u3d = .false.
+ dta_bdy(ib_bdy)%ll_v3d = .false.
+ CASE('zerograd')
+ IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities'
+ dta_bdy(ib_bdy)%ll_u3d = .false.
+ dta_bdy(ib_bdy)%ll_v3d = .false.
+ CASE('zero')
+ IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)'
+ dta_bdy(ib_bdy)%ll_u3d = .false.
+ dta_bdy(ib_bdy)%ll_v3d = .false.
+ CASE('orlanski')
+ IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_u3d = .true.
+ dta_bdy(ib_bdy)%ll_v3d = .true.
+ CASE('orlanski_npo')
+ IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_u3d = .true.
+ dta_bdy(ib_bdy)%ll_v3d = .true.
+ CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
+ END SELECT
+ IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+ SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !
+ CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
+ CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
+ CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
+ END SELECT
+ ENDIF
+
+ IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
+ IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
+ IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
+ ln_dyn3d_dmp(ib_bdy)=.false.
+ ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
+ CALL ctl_stop( 'Use FRS OR relaxation' )
+ ELSE
+ IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone'
+ IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+ IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+ dta_bdy(ib_bdy)%ll_u3d = .true.
+ dta_bdy(ib_bdy)%ll_v3d = .true.
+ ENDIF
+ ELSE
+ IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'
+ ENDIF
+ IF(lwp) WRITE(numout,*)
+
+ IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
+ SELECT CASE( cn_tra(ib_bdy) )
+ CASE('none')
+ IF(lwp) WRITE(numout,*) ' no open boundary condition'
+ dta_bdy(ib_bdy)%ll_tem = .false.
+ dta_bdy(ib_bdy)%ll_sal = .false.
+ CASE('frs')
+ IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
+ dta_bdy(ib_bdy)%ll_tem = .true.
+ dta_bdy(ib_bdy)%ll_sal = .true.
+ CASE('specified')
+ IF(lwp) WRITE(numout,*) ' Specified value'
+ dta_bdy(ib_bdy)%ll_tem = .true.
+ dta_bdy(ib_bdy)%ll_sal = .true.
+ CASE('neumann')
+ IF(lwp) WRITE(numout,*) ' Neumann conditions'
+ dta_bdy(ib_bdy)%ll_tem = .false.
+ dta_bdy(ib_bdy)%ll_sal = .false.
+ CASE('runoff')
+ IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity'
+ dta_bdy(ib_bdy)%ll_tem = .false.
+ dta_bdy(ib_bdy)%ll_sal = .false.
+ CASE('orlanski')
+ IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_tem = .true.
+ dta_bdy(ib_bdy)%ll_sal = .true.
+ CASE('orlanski_npo')
+ IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
+ dta_bdy(ib_bdy)%ll_tem = .true.
+ dta_bdy(ib_bdy)%ll_sal = .true.
+ CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' )
+ END SELECT
+ IF( cn_tra(ib_bdy) /= 'none' ) THEN
+ SELECT CASE( nn_tra_dta(ib_bdy) ) !
+ CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
+ CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
+ CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
+ END SELECT
+ ENDIF
+
+ IF ( ln_tra_dmp(ib_bdy) ) THEN
+ IF ( cn_tra(ib_bdy) == 'none' ) THEN
+ IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
+ ln_tra_dmp(ib_bdy)=.false.
+ ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
+ CALL ctl_stop( 'Use FRS OR relaxation' )
+ ELSE
+ IF(lwp) WRITE(numout,*) ' + T/S relaxation zone'
+ IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+ IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
+ IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+ dta_bdy(ib_bdy)%ll_tem = .true.
+ dta_bdy(ib_bdy)%ll_sal = .true.
+ ENDIF
+ ELSE
+ IF(lwp) WRITE(numout,*) ' NO T/S relaxation'
+ ENDIF
+ IF(lwp) WRITE(numout,*)
+
+#if defined key_si3
+ IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
+ SELECT CASE( cn_ice(ib_bdy) )
+ CASE('none')
+ IF(lwp) WRITE(numout,*) ' no open boundary condition'
+ dta_bdy(ib_bdy)%ll_a_i = .false.
+ dta_bdy(ib_bdy)%ll_h_i = .false.
+ dta_bdy(ib_bdy)%ll_h_s = .false.
+ CASE('frs')
+ IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
+ dta_bdy(ib_bdy)%ll_a_i = .true.
+ dta_bdy(ib_bdy)%ll_h_i = .true.
+ dta_bdy(ib_bdy)%ll_h_s = .true.
+ CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' )
+ END SELECT
+ IF( cn_ice(ib_bdy) /= 'none' ) THEN
+ SELECT CASE( nn_ice_dta(ib_bdy) ) !
+ CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
+ CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
+ CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
+ END SELECT
+ ENDIF
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)
+ IF(lwp) WRITE(numout,*) ' sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)
+ IF(lwp) WRITE(numout,*) ' age of bdy sea-ice = ', rn_ice_age(ib_bdy)
+#endif
+
+ IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)
+ IF(lwp) WRITE(numout,*)
+ !
+ END DO
+
+ IF( nb_bdy > 0 ) THEN
+ IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)
+ IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
+ IF(lwp) WRITE(numout,*)
+ SELECT CASE ( nn_volctl )
+ CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'
+ CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'
+ CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )
+ END SELECT
+ IF(lwp) WRITE(numout,*)
+ ELSE
+ IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
+ IF(lwp) WRITE(numout,*)
+ ENDIF
+ IF( nb_jpk_bdy > 0 ) THEN
+ IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
+ ELSE
+ IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'
+ ENDIF
+ ENDIF
+
+ ! -------------------------------------------------
+ ! Initialise indices arrays for open boundaries
+ ! -------------------------------------------------
+
+ ! Work out global dimensions of boundary data
+ ! ---------------------------------------------
+ REWIND( numnam_cfg )
+
+ nblendta(:,:) = 0
+ nbdysege = 0
+ nbdysegw = 0
+ nbdysegn = 0
+ nbdysegs = 0
+ icount = 0 ! count user defined segments
+ ! Dimensions below are used to allocate arrays to read external data
+ jpbdtas = 1 ! Maximum size of boundary data (structured case)
+ jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
+
+ DO ib_bdy = 1, nb_bdy
+
+ IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
+
+ icount = icount + 1
+ ! No REWIND here because may need to read more than one nambdy_index namelist.
+ ! Read only namelist_cfg to avoid unseccessfull overwrite
+ ! keep full control of the configuration namelist
+ READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
+904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, nambdy_index )
+
+ SELECT CASE ( TRIM(ctypebdy) )
+ CASE( 'N' )
+ IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
+ nbdyind = jpjglo - 2 ! set boundary to whole side of model domain.
+ nbdybeg = 2
+ nbdyend = jpiglo - 1
+ ENDIF
+ nbdysegn = nbdysegn + 1
+ npckgn(nbdysegn) = ib_bdy ! Save bdy package number
+ jpjnob(nbdysegn) = nbdyind
+ jpindt(nbdysegn) = nbdybeg
+ jpinft(nbdysegn) = nbdyend
+ !
+ CASE( 'S' )
+ IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
+ nbdyind = 2 ! set boundary to whole side of model domain.
+ nbdybeg = 2
+ nbdyend = jpiglo - 1
+ ENDIF
+ nbdysegs = nbdysegs + 1
+ npckgs(nbdysegs) = ib_bdy ! Save bdy package number
+ jpjsob(nbdysegs) = nbdyind
+ jpisdt(nbdysegs) = nbdybeg
+ jpisft(nbdysegs) = nbdyend
+ !
+ CASE( 'E' )
+ IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
+ nbdyind = jpiglo - 2 ! set boundary to whole side of model domain.
+ nbdybeg = 2
+ nbdyend = jpjglo - 1
+ ENDIF
+ nbdysege = nbdysege + 1
+ npckge(nbdysege) = ib_bdy ! Save bdy package number
+ jpieob(nbdysege) = nbdyind
+ jpjedt(nbdysege) = nbdybeg
+ jpjeft(nbdysege) = nbdyend
+ !
+ CASE( 'W' )
+ IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
+ nbdyind = 2 ! set boundary to whole side of model domain.
+ nbdybeg = 2
+ nbdyend = jpjglo - 1
+ ENDIF
+ nbdysegw = nbdysegw + 1
+ npckgw(nbdysegw) = ib_bdy ! Save bdy package number
+ jpiwob(nbdysegw) = nbdyind
+ jpjwdt(nbdysegw) = nbdybeg
+ jpjwft(nbdysegw) = nbdyend
+ !
+ CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
+ END SELECT
+
+ ! For simplicity we assume that in case of straight bdy, arrays have the same length
+ ! (even if it is true that last tangential velocity points
+ ! are useless). This simplifies a little bit boundary data format (and agrees with format
+ ! used so far in obc package)
+
+ nblendta(1:jpbgrd,ib_bdy) = (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
+ jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
+ IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
+ & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
+
+ ELSE ! Read size of arrays in boundary coordinates file.
+ CALL iom_open( cn_coords_file(ib_bdy), inum )
+ DO igrd = 1, jpbgrd
+ id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )
+ nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
+ jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
+ END DO
+ CALL iom_close( inum )
+ !
+ ENDIF
+ !
+ END DO ! ib_bdy
+
+ IF (nb_bdy>0) THEN
+ jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
+
+ ! Allocate arrays
+ !---------------
+ ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
+ & nbrdta(jpbdta, jpbgrd, nb_bdy) )
+
+ IF( nb_jpk_bdy>0 ) THEN
+ ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
+ ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
+ ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
+ ELSE
+ ALLOCATE( dta_global(jpbdtau, 1, jpk) )
+ ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
+ ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
+ ENDIF
+
+ IF ( icount>0 ) THEN
+ IF( nb_jpk_bdy>0 ) THEN
+ ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
+ ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
+ ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
+ ELSE
+ ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
+ ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
+ ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO
+ ENDIF
+ ENDIF
+ !
+ ENDIF
+
+ ! Now look for crossings in user (namelist) defined open boundary segments:
+ !--------------------------------------------------------------------------
+ IF( icount>0 ) CALL bdy_ctl_seg
+
+ ! Calculate global boundary index arrays or read in from file
+ !------------------------------------------------------------
+ ! 1. Read global index arrays from boundary coordinates file.
+ DO ib_bdy = 1, nb_bdy
+ !
+ IF( ln_coords_file(ib_bdy) ) THEN
+ !
+ CALL iom_open( cn_coords_file(ib_bdy), inum )
+ DO igrd = 1, jpbgrd
+ CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+ DO ii = 1,nblendta(igrd,ib_bdy)
+ nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+ END DO
+ CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+ DO ii = 1,nblendta(igrd,ib_bdy)
+ nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+ END DO
+ CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+ DO ii = 1,nblendta(igrd,ib_bdy)
+ nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+ END DO
+ !
+ ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
+ IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
+ IF (ibr_max < nn_rimwidth(ib_bdy)) &
+ CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
+ END DO
+ CALL iom_close( inum )
+ !
+ ENDIF
+ !
+ END DO
+
+ ! 2. Now fill indices corresponding to straight open boundary arrays:
+ ! East
+ !-----
+ DO iseg = 1, nbdysege
+ ib_bdy = npckge(iseg)
+ !
+ ! ------------ T points -------------
+ igrd=1
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ij = jpjedt(iseg), jpjeft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ U points -------------
+ igrd=2
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ij = jpjedt(iseg), jpjeft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ V points -------------
+ igrd=3
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+! DO ij = jpjedt(iseg), jpjeft(iseg) - 1
+ DO ij = jpjedt(iseg), jpjeft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ ENDDO
+ ENDDO
+ !
+ ! West
+ !-----
+ DO iseg = 1, nbdysegw
+ ib_bdy = npckgw(iseg)
+ !
+ ! ------------ T points -------------
+ igrd=1
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ij = jpjwdt(iseg), jpjwft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ U points -------------
+ igrd=2
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ij = jpjwdt(iseg), jpjwft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ V points -------------
+ igrd=3
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
+ DO ij = jpjwdt(iseg), jpjwft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+ nbjdta(icount, igrd, ib_bdy) = ij
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ ENDDO
+ ENDDO
+ !
+ ! North
+ !-----
+ DO iseg = 1, nbdysegn
+ ib_bdy = npckgn(iseg)
+ !
+ ! ------------ T points -------------
+ igrd=1
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ii = jpindt(iseg), jpinft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ U points -------------
+ igrd=2
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+! DO ii = jpindt(iseg), jpinft(iseg) - 1
+ DO ii = jpindt(iseg), jpinft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ ENDDO
+ !
+ ! ------------ V points -------------
+ igrd=3
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ii = jpindt(iseg), jpinft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ ENDDO
+ !
+ ! South
+ !-----
+ DO iseg = 1, nbdysegs
+ ib_bdy = npckgs(iseg)
+ !
+ ! ------------ T points -------------
+ igrd=1
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ii = jpisdt(iseg), jpisft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ !
+ ! ------------ U points -------------
+ igrd=2
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+! DO ii = jpisdt(iseg), jpisft(iseg) - 1
+ DO ii = jpisdt(iseg), jpisft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+ ENDDO
+ !
+ ! ------------ V points -------------
+ igrd=3
+ icount=0
+ DO ir = 1, nn_rimwidth(ib_bdy)
+ DO ii = jpisdt(iseg), jpisft(iseg)
+ icount = icount + 1
+ nbidta(icount, igrd, ib_bdy) = ii
+ nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+ nbrdta(icount, igrd, ib_bdy) = ir
+ ENDDO
+ ENDDO
+ ENDDO
+
+ ! Deal with duplicated points
+ !-----------------------------
+ ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
+ ! if their distance to the bdy is greater than the other
+ ! If their distance are the same, just keep only one to avoid updating a point twice
+ DO igrd = 1, jpbgrd
+ DO ib_bdy1 = 1, nb_bdy
+ DO ib_bdy2 = 1, nb_bdy
+ IF (ib_bdy1/=ib_bdy2) THEN
+ DO ib1 = 1, nblendta(igrd,ib_bdy1)
+ DO ib2 = 1, nblendta(igrd,ib_bdy2)
+ IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
+ & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
+! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
+! & nbidta(ib1, igrd, ib_bdy1), &
+! & nbjdta(ib2, igrd, ib_bdy2)
+ ! keep only points with the lowest distance to boundary:
+ IF (nbrdta(ib1, igrd, ib_bdy1)nbrdta(ib2, igrd, ib_bdy2)) THEN
+ nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+ nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+ ! Arbitrary choice if distances are the same:
+ ELSE
+ nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+ nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+ ENDIF
+ END IF
+ END DO
+ END DO
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ ! Work out dimensions of boundary data on each processor
+ ! ------------------------------------------------------
+
+ ! Rather assume that boundary data indices are given on global domain
+ ! TO BE DISCUSSED ?
+! iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
+! ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
+! is = mjg(1) + 1 ! if monotasking and no zoom, is=2
+! in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
+ iwe = mig(1) - 1 + 2 ! if monotasking and no zoom, iw=2
+ ies = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
+ iso = mjg(1) - 1 + 2 ! if monotasking and no zoom, is=2
+ ino = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
+
+ ALLOCATE( nbondi_bdy(nb_bdy))
+ ALLOCATE( nbondj_bdy(nb_bdy))
+ nbondi_bdy(:)=2
+ nbondj_bdy(:)=2
+ ALLOCATE( nbondi_bdy_b(nb_bdy))
+ ALLOCATE( nbondj_bdy_b(nb_bdy))
+ nbondi_bdy_b(:)=2
+ nbondj_bdy_b(:)=2
+
+ ! Work out dimensions of boundary data on each neighbour process
+ IF(nbondi == 0) THEN
+ iw_b(1) = 1 + nimppt(nowe+1)
+ ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
+ is_b(1) = 1 + njmppt(nowe+1)
+ in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
+
+ iw_b(2) = 1 + nimppt(noea+1)
+ ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
+ is_b(2) = 1 + njmppt(noea+1)
+ in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
+ ELSEIF(nbondi == 1) THEN
+ iw_b(1) = 1 + nimppt(nowe+1)
+ ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
+ is_b(1) = 1 + njmppt(nowe+1)
+ in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
+ ELSEIF(nbondi == -1) THEN
+ iw_b(2) = 1 + nimppt(noea+1)
+ ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
+ is_b(2) = 1 + njmppt(noea+1)
+ in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
+ ENDIF
+
+ IF(nbondj == 0) THEN
+ iw_b(3) = 1 + nimppt(noso+1)
+ ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
+ is_b(3) = 1 + njmppt(noso+1)
+ in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
+
+ iw_b(4) = 1 + nimppt(nono+1)
+ ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
+ is_b(4) = 1 + njmppt(nono+1)
+ in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
+ ELSEIF(nbondj == 1) THEN
+ iw_b(3) = 1 + nimppt(noso+1)
+ ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
+ is_b(3) = 1 + njmppt(noso+1)
+ in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
+ ELSEIF(nbondj == -1) THEN
+ iw_b(4) = 1 + nimppt(nono+1)
+ ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
+ is_b(4) = 1 + njmppt(nono+1)
+ in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
+ ENDIF
+
+ DO ib_bdy = 1, nb_bdy
+ DO igrd = 1, jpbgrd
+ icount = 0
+ icountr = 0
+ idx_bdy(ib_bdy)%nblen(igrd) = 0
+ idx_bdy(ib_bdy)%nblenrim(igrd) = 0
+ DO ib = 1, nblendta(igrd,ib_bdy)
+ ! check that data is in correct order in file
+ ibm1 = MAX(1,ib-1)
+ IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
+ IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
+ CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
+ & ' in order of distance from edge nbr A utility for re-ordering ', &
+ & ' boundary coordinates and data files exists in the TOOLS/OBC directory')
+ ENDIF
+ ENDIF
+ ! check if point is in local domain
+ IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN
+ !
+ icount = icount + 1
+ !
+ IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
+ ENDIF
+ END DO
+ idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
+ idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
+ END DO ! igrd
+
+ ! Allocate index arrays for this boundary set
+ !--------------------------------------------
+ ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
+ ALLOCATE( idx_bdy(ib_bdy)%nbi (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbj (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbr (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) , &
+ & idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
+
+ ! Dispatch mapping indices and discrete distances on each processor
+ ! -----------------------------------------------------------------
+
+ com_east = 0
+ com_west = 0
+ com_south = 0
+ com_north = 0
+
+ com_east_b = 0
+ com_west_b = 0
+ com_south_b = 0
+ com_north_b = 0
+
+ DO igrd = 1, jpbgrd
+ icount = 0
+ ! Loop on rimwidth to ensure outermost points come first in the local arrays.
+ DO ir=1, nn_rimwidth(ib_bdy)
+ DO ib = 1, nblendta(igrd,ib_bdy)
+ ! check if point is in local domain and equals ir
+ IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ !
+ icount = icount + 1
+
+ ! Rather assume that boundary data indices are given on global domain
+ ! TO BE DISCUSSED ?
+! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
+! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
+ idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
+ idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
+ ! check if point has to be sent
+ ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
+ ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
+ if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then
+ com_east = 1
+ elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
+ com_west = 1
+ endif
+ if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
+ com_south = 1
+ elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then
+ com_north = 1
+ endif
+ idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
+ idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
+ ENDIF
+ ! check if point has to be received from a neighbour
+ IF(nbondi == 0) THEN
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
+ if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+ ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
+ if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
+ com_south = 1
+ elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
+ com_north = 1
+ endif
+ com_west_b = 1
+ endif
+ ENDIF
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
+ if((com_east_b .ne. 1) .and. (ii == 2)) then
+ ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
+ if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
+ com_south = 1
+ elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
+ com_north = 1
+ endif
+ com_east_b = 1
+ endif
+ ENDIF
+ ELSEIF(nbondi == 1) THEN
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
+ if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+ ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
+ if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
+ com_south = 1
+ elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
+ com_north = 1
+ endif
+ com_west_b = 1
+ endif
+ ENDIF
+ ELSEIF(nbondi == -1) THEN
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
+ if((com_east_b .ne. 1) .and. (ii == 2)) then
+ ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
+ if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
+ com_south = 1
+ elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
+ com_north = 1
+ endif
+ com_east_b = 1
+ endif
+ ENDIF
+ ENDIF
+ IF(nbondj == 0) THEN
+ IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
+ & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
+ com_north_b = 1
+ ENDIF
+ IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 &
+ &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
+ com_south_b = 1
+ ENDIF
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
+ if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
+ com_south_b = 1
+ endif
+ ENDIF
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
+ if((com_north_b .ne. 1) .and. (ij == 2)) then
+ com_north_b = 1
+ endif
+ ENDIF
+ ELSEIF(nbondj == 1) THEN
+ IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
+ & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
+ com_south_b = 1
+ ENDIF
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
+ if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
+ com_south_b = 1
+ endif
+ ENDIF
+ ELSEIF(nbondj == -1) THEN
+ IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
+ & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
+ com_north_b = 1
+ ENDIF
+ IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
+ & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
+ & nbrdta(ib,igrd,ib_bdy) == ir ) THEN
+ ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
+ if((com_north_b .ne. 1) .and. (ij == 2)) then
+ com_north_b = 1
+ endif
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+
+ ! definition of the i- and j- direction local boundaries arrays used for sending the boundaries
+ IF( (com_east == 1) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 0
+ ELSEIF( (com_east == 1) .and. (com_west == 0) ) THEN ; nbondi_bdy(ib_bdy) = -1
+ ELSEIF( (com_east == 0) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 1
+ ENDIF
+ IF( (com_north == 1) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 0
+ ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN ; nbondj_bdy(ib_bdy) = -1
+ ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 1
+ ENDIF
+
+ ! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries
+ IF( (com_east_b == 1) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 0
+ ELSEIF( (com_east_b == 1) .and. (com_west_b == 0) ) THEN ; nbondi_bdy_b(ib_bdy) = -1
+ ELSEIF( (com_east_b == 0) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 1
+ ENDIF
+ IF( (com_north_b == 1) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 0
+ ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN ; nbondj_bdy_b(ib_bdy) = -1
+ ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 1
+ ENDIF
+
+ ! Compute rim weights for FRS scheme
+ ! ----------------------------------
+ DO igrd = 1, jpbgrd
+ DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+ nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+ idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 ) ! tanh formulation
+! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
+! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)) ! linear
+ END DO
+ END DO
+
+ ! Compute damping coefficients
+ ! ----------------------------
+! DO igrd = 1, jpbgrd
+! DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+! nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+! idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
+! & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
+! idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
+! & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
+! END DO
+! END DO
+! Hard coded COMODO sponge:
+ DO igrd = 1, jpbgrd
+ DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+ nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+ idx_bdy(ib_bdy)%nbd(ib,igrd) = 1._wp / ( rn_time_dmp(ib_bdy) * rday ) &
+ & *EXP(-5.0E-5*1000._wp*(nbr-1._wp)) ! exponential
+ idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1._wp / ( rn_time_dmp_out(ib_bdy) * rday ) &
+ & *EXP(-5.0E-5*1000._wp*(nbr-1._wp)) ! exponential
+ END DO
+ END DO
+ END DO
+
+ ! ------------------------------------------------------
+ ! Initialise masks and find normal/tangential directions
+ ! ------------------------------------------------------
+
+ ! Read global 2D mask at T-points: bdytmask
+ ! -----------------------------------------
+ ! bdytmask = 1 on the computational domain AND on open boundaries
+ ! = 0 elsewhere
+
+ bdytmask(:,:) = ssmask(:,:)
+
+ ! Derive mask on U and V grid from mask on T grid
+
+ bdyumask(:,:) = 0._wp
+ bdyvmask(:,:) = 0._wp
+ DO ij = 1, jpjm1
+ DO ii = 1, jpim1
+ bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
+ bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1)
+ END DO
+ END DO
+ CALL lbc_lnk_multi( bdyumask, 'U', 1. , bdyvmask, 'V', 1. ) ! Lateral boundary cond.
+
+ ! bdy masks are now set to zero on boundary points:
+ !
+ igrd = 1
+ DO ib_bdy = 1, nb_bdy
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
+ END DO
+ END DO
+ !
+ igrd = 2
+ DO ib_bdy = 1, nb_bdy
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
+ END DO
+ END DO
+ !
+ igrd = 3
+ DO ib_bdy = 1, nb_bdy
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
+ END DO
+ END DO
+
+ ! For the flagu/flagv calculation below we require a version of fmask without
+ ! the land boundary condition (shlat) included:
+ zfmask(:,:) = 0
+ DO ij = 2, jpjm1
+ DO ii = 2, jpim1
+ zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) &
+ & * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
+ END DO
+ END DO
+
+ ! Lateral boundary conditions
+ CALL lbc_lnk( zfmask, 'F', 1. )
+ CALL lbc_lnk_multi( bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. )
+ DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
+
+ idx_bdy(ib_bdy)%flagu(:,:) = 0._wp
+ idx_bdy(ib_bdy)%flagv(:,:) = 0._wp
+ icount = 0
+
+ ! Calculate relationship of U direction to the local orientation of the boundary
+ ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
+ ! flagu = 0 : u is tangential
+ ! flagu = 1 : u is normal to the boundary and is direction is inward
+
+ DO igrd = 1, jpbgrd
+ SELECT CASE( igrd )
+ CASE( 1 ) ; pmask => umask (:,:,1) ; i_offset = 0
+ CASE( 2 ) ; pmask => bdytmask(:,:) ; i_offset = 1
+ CASE( 3 ) ; pmask => zfmask (:,:) ; i_offset = 0
+ END SELECT
+ icount = 0
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+ nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+ zefl = pmask(nbi+i_offset-1,nbj)
+ zwfl = pmask(nbi+i_offset,nbj)
+ ! This error check only works if you are using the bdyXmask arrays
+ IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN
+ icount = icount + 1
+ IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
+ ELSE
+ idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl
+ ENDIF
+ END DO
+ IF( icount /= 0 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
+ ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
+ IF(lwp) WRITE(numout,*) ' ========== '
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ END DO
+
+ ! Calculate relationship of V direction to the local orientation of the boundary
+ ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
+ ! flagv = 0 : v is tangential
+ ! flagv = 1 : v is normal to the boundary and is direction is inward
+
+ DO igrd = 1, jpbgrd
+ SELECT CASE( igrd )
+ CASE( 1 ) ; pmask => vmask (:,:,1) ; j_offset = 0
+ CASE( 2 ) ; pmask => zfmask(:,:) ; j_offset = 0
+ CASE( 3 ) ; pmask => bdytmask ; j_offset = 1
+ END SELECT
+ icount = 0
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+ nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+ znfl = pmask(nbi,nbj+j_offset-1)
+ zsfl = pmask(nbi,nbj+j_offset )
+ ! This error check only works if you are using the bdyXmask arrays
+ IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
+ IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
+ icount = icount + 1
+ ELSE
+ idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl
+ END IF
+ END DO
+ IF( icount /= 0 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
+ ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
+ IF(lwp) WRITE(numout,*) ' ========== '
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ END DO
+ !
+ END DO
+
+ ! Compute total lateral surface for volume correction:
+ ! ----------------------------------------------------
+ ! JC: this must be done at each time step with non-linear free surface
+ bdysurftot = 0._wp
+ IF( ln_vol ) THEN
+ igrd = 2 ! Lateral surface at U-points
+ DO ib_bdy = 1, nb_bdy
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+ nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+ flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
+ bdysurftot = bdysurftot + hu_n (nbi , nbj) &
+ & * e2u (nbi , nbj) * ABS( flagu ) &
+ & * tmask_i(nbi , nbj) &
+ & * tmask_i(nbi+1, nbj)
+ END DO
+ END DO
+
+ igrd=3 ! Add lateral surface at V-points
+ DO ib_bdy = 1, nb_bdy
+ DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+ nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+ flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
+ bdysurftot = bdysurftot + hv_n (nbi, nbj ) &
+ & * e1v (nbi, nbj ) * ABS( flagv ) &
+ & * tmask_i(nbi, nbj ) &
+ & * tmask_i(nbi, nbj+1)
+ END DO
+ END DO
+ !
+ IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
+ END IF
+ !
+ ! Tidy up
+ !--------
+ IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta )
+ !
+ END SUBROUTINE bdy_segs
+
+
+ SUBROUTINE bdy_ctl_seg
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE bdy_ctl_seg ***
+ !!
+ !! ** Purpose : Check straight open boundary segments location
+ !!
+ !! ** Method : - Look for open boundary corners
+ !! - Check that segments start or end on land
+ !!----------------------------------------------------------------------
+ INTEGER :: ib, ib1, ib2, ji ,jj, itest, icyclic
+ INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
+ REAL(wp), DIMENSION(2) :: ztestmask
+ !!----------------------------------------------------------------------
+ !
+ IF (lwp) WRITE(numout,*) ' '
+ IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
+ IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+ !
+ IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
+ IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
+ IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
+ IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
+ ! 1. Check bounds
+ !----------------
+ DO ib = 1, nbdysegn
+ IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
+ IF ((jpjnob(ib).ge.jpjglo-1).or.&
+ &(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
+ IF (jpindt(ib).gt.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+ IF (jpindt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
+ IF (jpinft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
+ END DO
+ !
+ DO ib = 1, nbdysegs
+ IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
+ IF ((jpjsob(ib).ge.jpjglo-1).or.&
+ &(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
+ IF (jpisdt(ib).gt.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+ IF (jpisdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
+ IF (jpisft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
+ END DO
+ !
+ DO ib = 1, nbdysege
+ IF (lwp) WRITE(numout,*) '**check east seg bounds pckg: ', npckge(ib)
+ IF ((jpieob(ib).ge.jpiglo-1).or.&
+ &(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
+ IF (jpjedt(ib).gt.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+ IF (jpjedt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
+ IF (jpjeft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
+ END DO
+ !
+ DO ib = 1, nbdysegw
+ IF (lwp) WRITE(numout,*) '**check west seg bounds pckg: ', npckgw(ib)
+ IF ((jpiwob(ib).ge.jpiglo-1).or.&
+ &(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
+ IF (jpjwdt(ib).gt.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+ IF (jpjwdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
+ IF (jpjwft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
+ ENDDO
+ !
+ !
+ ! 2. Look for segment crossings
+ !------------------------------
+ IF (lwp) WRITE(numout,*) '**Look for segments corners :'
+ !
+ itest = 0 ! corner number
+ !
+ ! flag to detect if start or end of open boundary belongs to a corner
+ ! if not (=0), it must be on land.
+ ! if a corner is detected, save bdy package number for further tests
+ icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
+ ! South/West crossings
+ IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
+ DO ib1 = 1, nbdysegw
+ DO ib2 = 1, nbdysegs
+ IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
+ & ( jpisft(ib2)>=jpiwob(ib1)).AND. &
+ & ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
+ & ( jpjwft(ib1)>=jpjsob(ib2))) THEN
+ IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
+ ! We have a possible South-West corner
+! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
+! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
+ icornw(ib1,1) = npckgs(ib2)
+ icorns(ib2,1) = npckgw(ib1)
+ ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+ & jpisft(ib2), jpjwft(ib1)
+ IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
+ IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
+ & ' and South segment: ',npckgs(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
+ IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1) , &
+ & ' and South segment: ',npckgs(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop+1
+ END IF
+ END IF
+ END DO
+ END DO
+ END IF
+ !
+ ! South/East crossings
+ IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
+ DO ib1 = 1, nbdysege
+ DO ib2 = 1, nbdysegs
+ IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
+ & ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
+ & ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
+ & ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
+ IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
+ ! We have a possible South-East corner
+! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
+! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
+ icorne(ib1,1) = npckgs(ib2)
+ icorns(ib2,2) = npckge(ib1)
+ ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+ & jpisdt(ib2), jpjeft(ib1)
+ IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
+ IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
+ & ' and South segment: ',npckgs(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
+ IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
+ & ' and South segment: ',npckgs(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ END IF
+ END IF
+ END DO
+ END DO
+ END IF
+ !
+ ! North/West crossings
+ IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
+ DO ib1 = 1, nbdysegw
+ DO ib2 = 1, nbdysegn
+ IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
+ & ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
+ & ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
+ & ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
+ IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
+ ! We have a possible North-West corner
+! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
+! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
+ icornw(ib1,2) = npckgn(ib2)
+ icornn(ib2,1) = npckgw(ib1)
+ ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+ & jpinft(ib2), jpjwdt(ib1)
+ IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
+ IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
+ & ' and North segment: ',npckgn(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
+ IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1), &
+ & ' and North segment: ',npckgn(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ END IF
+ END IF
+ END DO
+ END DO
+ END IF
+ !
+ ! North/East crossings
+ IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
+ DO ib1 = 1, nbdysege
+ DO ib2 = 1, nbdysegn
+ IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
+ & ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
+ & ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
+ & ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
+ IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
+ ! We have a possible North-East corner
+! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
+! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
+ icorne(ib1,2) = npckgn(ib2)
+ icornn(ib2,2) = npckge(ib1)
+ ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+ & jpindt(ib2), jpjedt(ib1)
+ IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
+ IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
+ & ' and North segment: ',npckgn(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
+ IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
+ & ' and North segment: ',npckgn(ib2)
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ END IF
+ END IF
+ END DO
+ END DO
+ END IF
+ !
+ ! 3. Check if segment extremities are on land
+ !--------------------------------------------
+ !
+ icyclic=0
+ ! West segments
+ DO ib = 1, nbdysegw
+ ! get mask at boundary extremities:
+ ztestmask(1:2)=0.
+ DO ji = 1, jpi
+ DO jj = 1, jpj
+ IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
+ & ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
+ IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
+ & ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
+
+ IF (ztestmask(1)==1) THEN
+ IF (icornw(ib,1)==0) THEN
+ IF ((jperio==2).AND.(jpjwdt(ib)==2)) THEN
+ ! Open boundary segment can be Cyclic North-South
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic N-S
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ELSE
+ ! This is a corner
+ IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
+ CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
+ itest=itest+1
+ ENDIF
+ ENDIF
+ IF (ztestmask(2)==1) THEN
+ IF (icornw(ib,2)==0) THEN
+ IF ((jperio==2).AND.(jpjwft(ib)==jpjglo-1)) THEN
+ ! Open boundary segment can be Cyclic North-South
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic N-S
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ELSE
+ ! This is a corner
+ IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
+ CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
+ itest=itest+1
+ ENDIF
+ ENDIF
+ IF ( icyclic==2 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' West open boundary segment ', npckgw(ib)
+ IF(lwp) WRITE(numout,*) ' is cyclic North-South'
+ IF(lwp) WRITE(numout,*)
+ ENDIF
+ END DO
+ !
+ icyclic=0
+ ! East segments
+ DO ib = 1, nbdysege
+ ! get mask at boundary extremities:
+ ztestmask(1:2)=0.
+ DO ji = 1, jpi
+ DO jj = 1, jpj
+ IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
+ & ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
+ IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
+ & ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
+
+ IF (ztestmask(1)==1) THEN
+ IF (icorne(ib,1)==0) THEN
+ IF ((jperio==2).AND.(jpjedt(ib)==2)) THEN
+ ! Open boundary segment can be Cyclic North-South
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic N-S
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ELSE
+ ! This is a corner
+ IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
+ CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
+ itest=itest+1
+ ENDIF
+ ENDIF
+ IF (ztestmask(2)==1) THEN
+ IF (icorne(ib,2)==0) THEN
+ IF ((jperio==2).AND.(jpjeft(ib)==jpjglo-1)) THEN
+ ! Open boundary segment can be Cyclic North-South
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic N-S
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ELSE
+ ! This is a corner
+ IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
+ CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
+ itest=itest+1
+ ENDIF
+ ENDIF
+ IF ( icyclic==2 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' East open boundary segment ', npckge(ib)
+ IF(lwp) WRITE(numout,*) ' is cyclic North-South'
+ IF(lwp) WRITE(numout,*)
+ ENDIF
+ END DO
+ !
+ ! South segments
+ icyclic=0
+ DO ib = 1, nbdysegs
+ ! get mask at boundary extremities:
+ ztestmask(1:2)=0.
+ DO ji = 1, jpi
+ DO jj = 1, jpj
+ IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
+ & ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
+ IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
+ & ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
+
+ IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
+ IF ((jperio==1).AND.(jpisdt(ib)==2)) THEN
+ ! Open boundary segment can be Cyclic East-West
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic east-west
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ENDIF
+ IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
+ IF ((jperio==1).AND.(jpisft(ib)==jpiglo-1)) THEN
+ ! Open boundary segment can be Cyclic East-West
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic east-west
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ENDIF
+ IF ( icyclic==2 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' South open boundary segment ', npckgs(ib)
+ IF(lwp) WRITE(numout,*) ' is cyclic East-West'
+ IF(lwp) WRITE(numout,*)
+ ENDIF
+ END DO
+ !
+ icyclic=0
+ ! North segments
+ DO ib = 1, nbdysegn
+ ! get mask at boundary extremities:
+ ztestmask(1:2)=0.
+ DO ji = 1, jpi
+ DO jj = 1, jpj
+ IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
+ & ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
+ IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
+ & ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)
+ END DO
+ END DO
+ IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
+
+ IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
+ IF ((jperio==1).AND.(jpindt(ib)==2)) THEN
+ ! Open boundary segment can be Cyclic East-West
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic east-west
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not start on land'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ENDIF
+ IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
+ IF ((jperio==1).AND.(jpindt(ib)==2)) THEN
+ ! Open boundary segment can be Cyclic East-West
+ icyclic = icyclic+1 ! Start indice is fine for a cyclic east-west
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
+ IF(lwp) WRITE(numout,*) ' ========== does not end on land'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ ENDIF
+ IF ( icyclic==2 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' North open boundary segment ', npckgn(ib)
+ IF(lwp) WRITE(numout,*) ' is cyclic East-West'
+ IF(lwp) WRITE(numout,*)
+ ENDIF
+ END DO
+ !
+ IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
+ !
+ ! Other tests TBD:
+ ! segments completly on land
+ ! optimized open boundary array length according to landmask
+ ! Nudging layers that overlap with interior domain
+ !
+ END SUBROUTINE bdy_ctl_seg
+
+
+ SUBROUTINE bdy_ctl_corn( ib1, ib2 )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE bdy_ctl_corn ***
+ !!
+ !! ** Purpose : Check numerical schemes consistency between
+ !! segments having a common corner
+ !!
+ !! ** Method :
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT(in) :: ib1, ib2
+ INTEGER :: itest
+ !!----------------------------------------------------------------------
+ itest = 0
+
+ IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
+ IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
+ IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
+ !
+ IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
+ IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
+ IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
+ !
+ IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
+ !
+ IF( itest>0 ) THEN
+ IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
+ IF(lwp) WRITE(numout,*) ' ========== have different open bdy schemes'
+ IF(lwp) WRITE(numout,*)
+ nstop = nstop + 1
+ ENDIF
+ !
+ END SUBROUTINE bdy_ctl_corn
+
+ !!=================================================================================
+END MODULE bdyini
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/eosbn2.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/eosbn2.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/eosbn2.F90 (revision 10118)
@@ -0,0 +1,1705 @@
+MODULE eosbn2
+ !!==============================================================================
+ !! *** MODULE eosbn2 ***
+ !! Equation Of Seawater : in situ density - Brunt-Vaisala frequency
+ !!==============================================================================
+ !! History : OPA ! 1989-03 (O. Marti) Original code
+ !! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
+ !! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
+ !! 7.0 ! 1996-01 (G. Madec) statement function for e3
+ !! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
+ !! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
+ !! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
+ !! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
+ !! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
+ !! - ! 2003-08 (G. Madec) F90, free form
+ !! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function)
+ !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
+ !! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp
+ !! 3.7 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation
+ !! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state
+ !! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module
+ !! - ! 2014-09 (F. Roquet) add TEOS-10, S-EOS, and modify EOS-80
+ !! - ! 2015-06 (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! eos : generic interface of the equation of state
+ !! eos_insitu : Compute the in situ density
+ !! eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass
+ !! eos_insitu_2d : Compute the in situ density for 2d fields
+ !! bn2 : Compute the Brunt-Vaisala frequency
+ !! eos_rab : generic interface of in situ thermal/haline expansion ratio
+ !! eos_rab_3d : compute in situ thermal/haline expansion ratio
+ !! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields
+ !! eos_fzp_2d : freezing temperature for 2d fields
+ !! eos_fzp_0d : freezing temperature for scalar
+ !! eos_init : set eos parameters (namelist)
+ !!----------------------------------------------------------------------
+ USE dom_oce ! ocean space and time domain
+ USE phycst ! physical constants
+ USE stopar ! Stochastic T/S fluctuations
+ USE stopts ! Stochastic T/S fluctuations
+ !
+ USE in_out_manager ! I/O manager
+ USE lib_mpp ! MPP library
+ USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+ USE prtctl ! Print control
+ USE lbclnk ! ocean lateral boundary conditions
+ USE timing ! Timing
+
+ IMPLICIT NONE
+ PRIVATE
+
+ ! !! * Interface
+ INTERFACE eos
+ MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
+ END INTERFACE
+ !
+ INTERFACE eos_rab
+ MODULE PROCEDURE rab_3d, rab_2d, rab_0d
+ END INTERFACE
+ !
+ INTERFACE eos_fzp
+ MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
+ END INTERFACE
+ !
+ PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules
+ PUBLIC bn2 ! called by step module
+ PUBLIC eos_rab ! called by ldfslp, zdfddm, trabbl
+ PUBLIC eos_pt_from_ct ! called by sbcssm
+ PUBLIC eos_fzp ! called by traadv_cen2 and sbcice_... modules
+ PUBLIC eos_pen ! used for pe diagnostics in trdpen module
+ PUBLIC eos_init ! called by istate module
+
+ ! !!** Namelist nameos **
+ LOGICAL , PUBLIC :: ln_TEOS10 ! determine if eos_pt_from_ct is used to compute sst_m
+ LOGICAL , PUBLIC :: ln_EOS80 ! determine if eos_pt_from_ct is used to compute sst_m
+ LOGICAL , PUBLIC :: ln_SEOS ! determine if eos_pt_from_ct is used to compute sst_m
+
+ ! Parameters
+ LOGICAL , PUBLIC :: l_useCT ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise
+ INTEGER , PUBLIC :: neos ! Identifier for equation of state used
+
+ INTEGER , PARAMETER :: np_teos10 = -1 ! parameter for using TEOS10
+ INTEGER , PARAMETER :: np_eos80 = 0 ! parameter for using EOS80
+ INTEGER , PARAMETER :: np_seos = 1 ! parameter for using Simplified Equation of state
+
+ ! !!! simplified eos coefficients (default value: Vallis 2006)
+ REAL(wp), PUBLIC :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff.
+ REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff.
+ REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2
+ REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2
+ REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T
+ REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S
+ REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt
+
+ ! TEOS10/EOS80 parameters
+ REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS
+
+ ! EOS parameters
+ REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
+ REAL(wp) :: EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510
+ REAL(wp) :: EOS020 , EOS120 , EOS220 , EOS320 , EOS420
+ REAL(wp) :: EOS030 , EOS130 , EOS230 , EOS330
+ REAL(wp) :: EOS040 , EOS140 , EOS240
+ REAL(wp) :: EOS050 , EOS150
+ REAL(wp) :: EOS060
+ REAL(wp) :: EOS001 , EOS101 , EOS201 , EOS301 , EOS401
+ REAL(wp) :: EOS011 , EOS111 , EOS211 , EOS311
+ REAL(wp) :: EOS021 , EOS121 , EOS221
+ REAL(wp) :: EOS031 , EOS131
+ REAL(wp) :: EOS041
+ REAL(wp) :: EOS002 , EOS102 , EOS202
+ REAL(wp) :: EOS012 , EOS112
+ REAL(wp) :: EOS022
+ REAL(wp) :: EOS003 , EOS103
+ REAL(wp) :: EOS013
+
+ ! ALPHA parameters
+ REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
+ REAL(wp) :: ALP010 , ALP110 , ALP210 , ALP310 , ALP410
+ REAL(wp) :: ALP020 , ALP120 , ALP220 , ALP320
+ REAL(wp) :: ALP030 , ALP130 , ALP230
+ REAL(wp) :: ALP040 , ALP140
+ REAL(wp) :: ALP050
+ REAL(wp) :: ALP001 , ALP101 , ALP201 , ALP301
+ REAL(wp) :: ALP011 , ALP111 , ALP211
+ REAL(wp) :: ALP021 , ALP121
+ REAL(wp) :: ALP031
+ REAL(wp) :: ALP002 , ALP102
+ REAL(wp) :: ALP012
+ REAL(wp) :: ALP003
+
+ ! BETA parameters
+ REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500
+ REAL(wp) :: BET010 , BET110 , BET210 , BET310 , BET410
+ REAL(wp) :: BET020 , BET120 , BET220 , BET320
+ REAL(wp) :: BET030 , BET130 , BET230
+ REAL(wp) :: BET040 , BET140
+ REAL(wp) :: BET050
+ REAL(wp) :: BET001 , BET101 , BET201 , BET301
+ REAL(wp) :: BET011 , BET111 , BET211
+ REAL(wp) :: BET021 , BET121
+ REAL(wp) :: BET031
+ REAL(wp) :: BET002 , BET102
+ REAL(wp) :: BET012
+ REAL(wp) :: BET003
+
+ ! PEN parameters
+ REAL(wp) :: PEN000 , PEN100 , PEN200 , PEN300 , PEN400
+ REAL(wp) :: PEN010 , PEN110 , PEN210 , PEN310
+ REAL(wp) :: PEN020 , PEN120 , PEN220
+ REAL(wp) :: PEN030 , PEN130
+ REAL(wp) :: PEN040
+ REAL(wp) :: PEN001 , PEN101 , PEN201
+ REAL(wp) :: PEN011 , PEN111
+ REAL(wp) :: PEN021
+ REAL(wp) :: PEN002 , PEN102
+ REAL(wp) :: PEN012
+
+ ! ALPHA_PEN parameters
+ REAL(wp) :: APE000 , APE100 , APE200 , APE300
+ REAL(wp) :: APE010 , APE110 , APE210
+ REAL(wp) :: APE020 , APE120
+ REAL(wp) :: APE030
+ REAL(wp) :: APE001 , APE101
+ REAL(wp) :: APE011
+ REAL(wp) :: APE002
+
+ ! BETA_PEN parameters
+ REAL(wp) :: BPE000 , BPE100 , BPE200 , BPE300
+ REAL(wp) :: BPE010 , BPE110 , BPE210
+ REAL(wp) :: BPE020 , BPE120
+ REAL(wp) :: BPE030
+ REAL(wp) :: BPE001 , BPE101
+ REAL(wp) :: BPE011
+ REAL(wp) :: BPE002
+
+ !! * Substitutions
+# include "vectopt_loop_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id: eosbn2.F90 9757 2018-06-07 08:54:30Z smasson $
+ !! Software governed by the CeCILL licence (./LICENSE)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE eos_insitu( pts, prd, pdep )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_insitu ***
+ !!
+ !! ** Purpose : Compute the in situ density (ratio rho/rau0) from
+ !! potential temperature and salinity using an equation of state
+ !! selected in the nameos namelist
+ !!
+ !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
+ !! with prd in situ density anomaly no units
+ !! t TEOS10: CT or EOS80: PT Celsius
+ !! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu
+ !! z depth meters
+ !! rho in situ density kg/m^3
+ !! rau0 reference density kg/m^3
+ !!
+ !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
+ !! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
+ !!
+ !! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
+ !! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
+ !!
+ !! ln_seos : simplified equation of state
+ !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
+ !! linear case function of T only: rn_alpha<>0, other coefficients = 0
+ !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
+ !! Vallis like equation: use default values of coefficients
+ !!
+ !! ** Action : compute prd , the in situ density (no units)
+ !!
+ !! References : Roquet et al, Ocean Modelling, in preparation (2014)
+ !! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+ !! TEOS-10 Manual, 2010
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-]
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m]
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt , zh , zs , ztm ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('eos-insitu')
+ !
+ SELECT CASE( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ zh = pdep(ji,jj,jk) * r1_Z0 ! depth
+ zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ ztm = tmask(ji,jj,jk) ! tmask
+ !
+ zn3 = EOS013*zt &
+ & + EOS103*zs+EOS003
+ !
+ zn2 = (EOS022*zt &
+ & + EOS112*zs+EOS012)*zt &
+ & + (EOS202*zs+EOS102)*zs+EOS002
+ !
+ zn1 = (((EOS041*zt &
+ & + EOS131*zs+EOS031)*zt &
+ & + (EOS221*zs+EOS121)*zs+EOS021)*zt &
+ & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
+ & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+ !
+ zn0 = (((((EOS060*zt &
+ & + EOS150*zs+EOS050)*zt &
+ & + (EOS240*zs+EOS140)*zs+EOS040)*zt &
+ & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
+ & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
+ & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
+ & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked)
+ !
+ END DO
+ END DO
+ END DO
+ !
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ zt = pts (ji,jj,jk,jp_tem) - 10._wp
+ zs = pts (ji,jj,jk,jp_sal) - 35._wp
+ zh = pdep (ji,jj,jk)
+ ztm = tmask(ji,jj,jk)
+ !
+ zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
+ & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
+ & - rn_nu * zt * zs
+ !
+ prd(ji,jj,jk) = zn * r1_rau0 * ztm ! density anomaly (masked)
+ END DO
+ END DO
+ END DO
+ !
+ END SELECT
+ !
+ IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk )
+ !
+ IF( ln_timing ) CALL timing_stop('eos-insitu')
+ !
+ END SUBROUTINE eos_insitu
+
+
+ SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_insitu_pot ***
+ !!
+ !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the
+ !! potential volumic mass (Kg/m3) from potential temperature and
+ !! salinity fields using an equation of state selected in the
+ !! namelist.
+ !!
+ !! ** Action : - prd , the in situ density (no units)
+ !! - prhop, the potential volumic mass (Kg/m3)
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-]
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced)
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m]
+ !
+ INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
+ INTEGER :: jdof
+ REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('eos-pot')
+ !
+ SELECT CASE ( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ ! Stochastic equation of state
+ IF ( ln_sto_eos ) THEN
+ ALLOCATE(zn0_sto(1:2*nn_sto_eos))
+ ALLOCATE(zn_sto(1:2*nn_sto_eos))
+ ALLOCATE(zsign(1:2*nn_sto_eos))
+ DO jsmp = 1, 2*nn_sto_eos, 2
+ zsign(jsmp) = 1._wp
+ zsign(jsmp+1) = -1._wp
+ END DO
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ ! compute density (2*nn_sto_eos) times:
+ ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
+ ! (2) for t-dt, s-ds (with the opposite fluctuation)
+ DO jsmp = 1, nn_sto_eos*2
+ jdof = (jsmp + 1) / 2
+ zh = pdep(ji,jj,jk) * r1_Z0 ! depth
+ zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature
+ zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
+ zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity
+ ztm = tmask(ji,jj,jk) ! tmask
+ !
+ zn3 = EOS013*zt &
+ & + EOS103*zs+EOS003
+ !
+ zn2 = (EOS022*zt &
+ & + EOS112*zs+EOS012)*zt &
+ & + (EOS202*zs+EOS102)*zs+EOS002
+ !
+ zn1 = (((EOS041*zt &
+ & + EOS131*zs+EOS031)*zt &
+ & + (EOS221*zs+EOS121)*zs+EOS021)*zt &
+ & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
+ & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+ !
+ zn0_sto(jsmp) = (((((EOS060*zt &
+ & + EOS150*zs+EOS050)*zt &
+ & + (EOS240*zs+EOS140)*zs+EOS040)*zt &
+ & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
+ & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
+ & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
+ & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+ !
+ zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
+ END DO
+ !
+ ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
+ prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
+ DO jsmp = 1, nn_sto_eos*2
+ prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface
+ !
+ prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked)
+ END DO
+ prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
+ prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos
+ END DO
+ END DO
+ END DO
+ DEALLOCATE(zn0_sto,zn_sto,zsign)
+ ! Non-stochastic equation of state
+ ELSE
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ zh = pdep(ji,jj,jk) * r1_Z0 ! depth
+ zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ ztm = tmask(ji,jj,jk) ! tmask
+ !
+ zn3 = EOS013*zt &
+ & + EOS103*zs+EOS003
+ !
+ zn2 = (EOS022*zt &
+ & + EOS112*zs+EOS012)*zt &
+ & + (EOS202*zs+EOS102)*zs+EOS002
+ !
+ zn1 = (((EOS041*zt &
+ & + EOS131*zs+EOS031)*zt &
+ & + (EOS221*zs+EOS121)*zs+EOS021)*zt &
+ & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
+ & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+ !
+ zn0 = (((((EOS060*zt &
+ & + EOS150*zs+EOS050)*zt &
+ & + (EOS240*zs+EOS140)*zs+EOS040)*zt &
+ & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
+ & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
+ & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
+ & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface
+ !
+ prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked)
+ END DO
+ END DO
+ END DO
+ ENDIF
+
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ zt = pts (ji,jj,jk,jp_tem) - 10._wp
+ zs = pts (ji,jj,jk,jp_sal) - 35._wp
+ zh = pdep (ji,jj,jk)
+ ztm = tmask(ji,jj,jk)
+ ! ! potential density referenced at the surface
+ zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt &
+ & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs &
+ & - rn_nu * zt * zs
+ prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
+ ! ! density anomaly (masked)
+ zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
+ prd(ji,jj,jk) = zn * r1_rau0 * ztm
+ !
+ END DO
+ END DO
+ END DO
+ !
+ END SELECT
+ !
+ IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
+ !
+ IF( ln_timing ) CALL timing_stop('eos-pot')
+ !
+ END SUBROUTINE eos_insitu_pot
+
+
+ SUBROUTINE eos_insitu_2d( pts, pdep, prd )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_insitu_2d ***
+ !!
+ !! ** Purpose : Compute the in situ density (ratio rho/rau0) from
+ !! potential temperature and salinity using an equation of state
+ !! selected in the nameos namelist. * 2D field case
+ !!
+ !! ** Action : - prd , the in situ density (no units) (unmasked)
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
+ ! ! 2 : salinity [psu]
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt , zh , zs ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('eos2d')
+ !
+ prd(:,:) = 0._wp
+ !
+ SELECT CASE( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ DO jj = 1, jpjm1
+ DO ji = 1, fs_jpim1 ! vector opt.
+ !
+ zh = pdep(ji,jj) * r1_Z0 ! depth
+ zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ !
+ zn3 = EOS013*zt &
+ & + EOS103*zs+EOS003
+ !
+ zn2 = (EOS022*zt &
+ & + EOS112*zs+EOS012)*zt &
+ & + (EOS202*zs+EOS102)*zs+EOS002
+ !
+ zn1 = (((EOS041*zt &
+ & + EOS131*zs+EOS031)*zt &
+ & + (EOS221*zs+EOS121)*zs+EOS021)*zt &
+ & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
+ & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+ !
+ zn0 = (((((EOS060*zt &
+ & + EOS150*zs+EOS050)*zt &
+ & + (EOS240*zs+EOS140)*zs+EOS040)*zt &
+ & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
+ & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
+ & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
+ & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ prd(ji,jj) = zn * r1_rau0 - 1._wp ! unmasked in situ density anomaly
+ !
+ END DO
+ END DO
+ !
+ CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions
+ !
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ DO jj = 1, jpjm1
+ DO ji = 1, fs_jpim1 ! vector opt.
+ !
+ zt = pts (ji,jj,jp_tem) - 10._wp
+ zs = pts (ji,jj,jp_sal) - 35._wp
+ zh = pdep (ji,jj) ! depth at the partial step level
+ !
+ zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
+ & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
+ & - rn_nu * zt * zs
+ !
+ prd(ji,jj) = zn * r1_rau0 ! unmasked in situ density anomaly
+ !
+ END DO
+ END DO
+ !
+ CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions
+ !
+ END SELECT
+ !
+ IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
+ !
+ IF( ln_timing ) CALL timing_stop('eos2d')
+ !
+ END SUBROUTINE eos_insitu_2d
+
+
+ SUBROUTINE rab_3d( pts, pab )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE rab_3d ***
+ !!
+ !! ** Purpose : Calculates thermal/haline expansion ratio at T-points
+ !!
+ !! ** Method : calculates alpha / beta at T-points
+ !!
+ !! ** Action : - pab : thermal/haline expansion ratio at T-points
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt , zh , zs , ztm ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('rab_3d')
+ !
+ SELECT CASE ( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth
+ zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ ztm = tmask(ji,jj,jk) ! tmask
+ !
+ ! alpha
+ zn3 = ALP003
+ !
+ zn2 = ALP012*zt + ALP102*zs+ALP002
+ !
+ zn1 = ((ALP031*zt &
+ & + ALP121*zs+ALP021)*zt &
+ & + (ALP211*zs+ALP111)*zs+ALP011)*zt &
+ & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+ !
+ zn0 = ((((ALP050*zt &
+ & + ALP140*zs+ALP040)*zt &
+ & + (ALP230*zs+ALP130)*zs+ALP030)*zt &
+ & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
+ & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
+ & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
+ !
+ ! beta
+ zn3 = BET003
+ !
+ zn2 = BET012*zt + BET102*zs+BET002
+ !
+ zn1 = ((BET031*zt &
+ & + BET121*zs+BET021)*zt &
+ & + (BET211*zs+BET111)*zs+BET011)*zt &
+ & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+ !
+ zn0 = ((((BET050*zt &
+ & + BET140*zs+BET040)*zt &
+ & + (BET230*zs+BET130)*zs+BET030)*zt &
+ & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
+ & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
+ & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
+ !
+ END DO
+ END DO
+ END DO
+ !
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
+ zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
+ zh = gdept_n(ji,jj,jk) ! depth in meters at t-point
+ ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask
+ !
+ zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+ pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm ! alpha
+ !
+ zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+ pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm ! beta
+ !
+ END DO
+ END DO
+ END DO
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ IF(ln_ctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
+ & tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
+ !
+ IF( ln_timing ) CALL timing_stop('rab_3d')
+ !
+ END SUBROUTINE rab_3d
+
+
+ SUBROUTINE rab_2d( pts, pdep, pab )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE rab_2d ***
+ !!
+ !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
+ !!
+ !! ** Action : - pab : thermal/haline expansion ratio at T-points
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt , zh , zs ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('rab_2d')
+ !
+ pab(:,:,:) = 0._wp
+ !
+ SELECT CASE ( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ DO jj = 1, jpjm1
+ DO ji = 1, fs_jpim1 ! vector opt.
+ !
+ zh = pdep(ji,jj) * r1_Z0 ! depth
+ zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ !
+ ! alpha
+ zn3 = ALP003
+ !
+ zn2 = ALP012*zt + ALP102*zs+ALP002
+ !
+ zn1 = ((ALP031*zt &
+ & + ALP121*zs+ALP021)*zt &
+ & + (ALP211*zs+ALP111)*zs+ALP011)*zt &
+ & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+ !
+ zn0 = ((((ALP050*zt &
+ & + ALP140*zs+ALP040)*zt &
+ & + (ALP230*zs+ALP130)*zs+ALP030)*zt &
+ & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
+ & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
+ & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(ji,jj,jp_tem) = zn * r1_rau0
+ !
+ ! beta
+ zn3 = BET003
+ !
+ zn2 = BET012*zt + BET102*zs+BET002
+ !
+ zn1 = ((BET031*zt &
+ & + BET121*zs+BET021)*zt &
+ & + (BET211*zs+BET111)*zs+BET011)*zt &
+ & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+ !
+ zn0 = ((((BET050*zt &
+ & + BET140*zs+BET040)*zt &
+ & + (BET230*zs+BET130)*zs+BET030)*zt &
+ & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
+ & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
+ & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(ji,jj,jp_sal) = zn / zs * r1_rau0
+ !
+ !
+ END DO
+ END DO
+ ! ! Lateral boundary conditions
+ CALL lbc_lnk_multi( pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )
+ !
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ DO jj = 1, jpjm1
+ DO ji = 1, fs_jpim1 ! vector opt.
+ !
+ zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
+ zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
+ zh = pdep (ji,jj) ! depth at the partial step level
+ !
+ zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+ pab(ji,jj,jp_tem) = zn * r1_rau0 ! alpha
+ !
+ zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+ pab(ji,jj,jp_sal) = zn * r1_rau0 ! beta
+ !
+ END DO
+ END DO
+ ! ! Lateral boundary conditions
+ CALL lbc_lnk_multi( pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ IF(ln_ctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
+ & tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
+ !
+ IF( ln_timing ) CALL timing_stop('rab_2d')
+ !
+ END SUBROUTINE rab_2d
+
+
+ SUBROUTINE rab_0d( pts, pdep, pab )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE rab_0d ***
+ !!
+ !! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
+ !!
+ !! ** Action : - pab : thermal/haline expansion ratio at T-points
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), INTENT(in ) :: pdep ! depth [m]
+ REAL(wp), DIMENSION(jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio
+ !
+ REAL(wp) :: zt , zh , zs ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('rab_0d')
+ !
+ pab(:) = 0._wp
+ !
+ SELECT CASE ( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ !
+ zh = pdep * r1_Z0 ! depth
+ zt = pts (jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ !
+ ! alpha
+ zn3 = ALP003
+ !
+ zn2 = ALP012*zt + ALP102*zs+ALP002
+ !
+ zn1 = ((ALP031*zt &
+ & + ALP121*zs+ALP021)*zt &
+ & + (ALP211*zs+ALP111)*zs+ALP011)*zt &
+ & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+ !
+ zn0 = ((((ALP050*zt &
+ & + ALP140*zs+ALP040)*zt &
+ & + (ALP230*zs+ALP130)*zs+ALP030)*zt &
+ & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
+ & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
+ & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(jp_tem) = zn * r1_rau0
+ !
+ ! beta
+ zn3 = BET003
+ !
+ zn2 = BET012*zt + BET102*zs+BET002
+ !
+ zn1 = ((BET031*zt &
+ & + BET121*zs+BET021)*zt &
+ & + (BET211*zs+BET111)*zs+BET011)*zt &
+ & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+ !
+ zn0 = ((((BET050*zt &
+ & + BET140*zs+BET040)*zt &
+ & + (BET230*zs+BET130)*zs+BET030)*zt &
+ & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
+ & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
+ & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+ !
+ zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+ !
+ pab(jp_sal) = zn / zs * r1_rau0
+ !
+ !
+ !
+ CASE( np_seos ) !== simplified EOS ==!
+ !
+ zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
+ zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
+ zh = pdep ! depth at the partial step level
+ !
+ zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+ pab(jp_tem) = zn * r1_rau0 ! alpha
+ !
+ zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+ pab(jp_sal) = zn * r1_rau0 ! beta
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ IF( ln_timing ) CALL timing_stop('rab_0d')
+ !
+ END SUBROUTINE rab_0d
+
+
+ SUBROUTINE bn2( pts, pab, pn2 )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE bn2 ***
+ !!
+ !! ** Purpose : Compute the local Brunt-Vaisala frequency at the
+ !! time-step of the input arguments
+ !!
+ !! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
+ !! where alpha and beta are given in pab, and computed on T-points.
+ !! N.B. N^2 is set one for all to zero at jk=1 in istate module.
+ !!
+ !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1]
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2]
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zaw, zbw, zrw ! local scalars
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('bn2')
+ !
+ DO jk = 2, jpkm1 ! interior points only (2=< jk =< jpkm1 )
+ DO jj = 1, jpj ! surface and bottom value set to zero one for all in istate.F90
+ DO ji = 1, jpi
+ zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) &
+ & / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
+ !
+ zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw
+ zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
+ !
+ pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
+ & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) &
+ & / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
+ END DO
+ END DO
+ END DO
+ !
+ IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', kdim=jpk )
+ !
+ IF( ln_timing ) CALL timing_stop('bn2')
+ !
+ END SUBROUTINE bn2
+
+
+ FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_pt_from_ct ***
+ !!
+ !! ** Purpose : Compute pot.temp. from cons. temp. [Celsius]
+ !!
+ !! ** Method : rational approximation (5/3th order) of TEOS-10 algorithm
+ !! checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
+ !!
+ !! Reference : TEOS-10, UNESCO
+ !! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ctmp ! Cons. Temp [Celsius]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
+ ! Leave result array automatic rather than making explicitly allocated
+ REAL(wp), DIMENSION(jpi,jpj) :: ptmp ! potential temperature [Celsius]
+ !
+ INTEGER :: ji, jj ! dummy loop indices
+ REAL(wp) :: zt , zs , ztm ! local scalars
+ REAL(wp) :: zn , zd ! local scalars
+ REAL(wp) :: zdeltaS , z1_S0 , z1_T0
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('eos_pt_from_ct')
+ !
+ zdeltaS = 5._wp
+ z1_S0 = 0.875_wp/35.16504_wp
+ z1_T0 = 1._wp/40._wp
+ !
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ zt = ctmp (ji,jj) * z1_T0
+ zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
+ ztm = tmask(ji,jj,1)
+ !
+ zn = ((((-2.1385727895e-01_wp*zt &
+ & - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt &
+ & + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt &
+ & + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt &
+ & + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs &
+ & +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt &
+ & + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs &
+ & -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
+ !
+ zd = (2.0035003456_wp*zt &
+ & -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt &
+ & + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
+ !
+ ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
+ !
+ END DO
+ END DO
+ !
+ IF( ln_timing ) CALL timing_stop('eos_pt_from_ct')
+ !
+ END FUNCTION eos_pt_from_ct
+
+
+ SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_fzp ***
+ !!
+ !! ** Purpose : Compute the freezing point temperature [Celsius]
+ !!
+ !! ** Method : UNESCO freezing point (ptf) in Celsius is given by
+ !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
+ !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
+ !!
+ !! Reference : UNESCO tech. papers in the marine science no. 28. 1978
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(out ) :: ptf ! freezing temperature [Celsius]
+ !
+ INTEGER :: ji, jj ! dummy loop indices
+ REAL(wp) :: zt, zs, z1_S0 ! local scalars
+ !!----------------------------------------------------------------------
+ !
+ SELECT CASE ( neos )
+ !
+ CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
+ !
+ z1_S0 = 1._wp / 35.16504_wp
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity
+ ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+ & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
+ END DO
+ END DO
+ ptf(:,:) = ptf(:,:) * psal(:,:)
+ !
+ IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
+ !
+ CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
+ !
+ ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &
+ & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)
+ !
+ IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ END SUBROUTINE eos_fzp_2d
+
+
+ SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_fzp ***
+ !!
+ !! ** Purpose : Compute the freezing point temperature [Celsius]
+ !!
+ !! ** Method : UNESCO freezing point (ptf) in Celsius is given by
+ !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
+ !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
+ !!
+ !! Reference : UNESCO tech. papers in the marine science no. 28. 1978
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in ) :: psal ! salinity [psu]
+ REAL(wp), INTENT(in ), OPTIONAL :: pdep ! depth [m]
+ REAL(wp), INTENT(out) :: ptf ! freezing temperature [Celsius]
+ !
+ REAL(wp) :: zs ! local scalars
+ !!----------------------------------------------------------------------
+ !
+ SELECT CASE ( neos )
+ !
+ CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
+ !
+ zs = SQRT( ABS( psal ) / 35.16504_wp ) ! square root salinity
+ ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+ & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
+ ptf = ptf * psal
+ !
+ IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
+ !
+ CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
+ !
+ ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal ) &
+ & - 2.154996e-4_wp * psal ) * psal
+ !
+ IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ END SUBROUTINE eos_fzp_0d
+
+
+ SUBROUTINE eos_pen( pts, pab_pe, ppen )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_pen ***
+ !!
+ !! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
+ !!
+ !! ** Method : PE is defined analytically as the vertical
+ !! primitive of EOS times -g integrated between 0 and z>0.
+ !! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
+ !! = 1/z * /int_0^z rd dz - rd
+ !! where rd is the density anomaly (see eos_rhd function)
+ !! ab_pe are partial derivatives of PE anomaly with respect to T and S:
+ !! ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
+ !! ab_pe(2) = 1/(rau0 gz) * dPE/dS + drd/dS = d(pen)/dS
+ !!
+ !! ** Action : - pen : PE anomaly given at T-points
+ !! : - pab_pe : given at T-points
+ !! pab_pe(:,:,:,jp_tem) is alpha_pe
+ !! pab_pe(:,:,:,jp_sal) is beta_pe
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zt , zh , zs , ztm ! local scalars
+ REAL(wp) :: zn , zn0, zn1, zn2 ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('eos_pen')
+ !
+ SELECT CASE ( neos )
+ !
+ CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ !
+ zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth
+ zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
+ zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
+ ztm = tmask(ji,jj,jk) ! tmask
+ !
+ ! potential energy non-linear anomaly
+ zn2 = (PEN012)*zt &
+ & + PEN102*zs+PEN002
+ !
+ zn1 = ((PEN021)*zt &
+ & + PEN111*zs+PEN011)*zt &
+ & + (PEN201*zs+PEN101)*zs+PEN001
+ !
+ zn0 = ((((PEN040)*zt &
+ & + PEN130*zs+PEN030)*zt &
+ & + (PEN220*zs+PEN120)*zs+PEN020)*zt &
+ & + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt &
+ & + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
+ !
+ zn = ( zn2 * zh + zn1 ) * zh + zn0
+ !
+ ppen(ji,jj,jk) = zn * zh * r1_rau0 * ztm
+ !
+ ! alphaPE non-linear anomaly
+ zn2 = APE002
+ !
+ zn1 = (APE011)*zt &
+ & + APE101*zs+APE001
+ !
+ zn0 = (((APE030)*zt &
+ & + APE120*zs+APE020)*zt &
+ & + (APE210*zs+APE110)*zs+APE010)*zt &
+ & + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
+ !
+ zn = ( zn2 * zh + zn1 ) * zh + zn0
+ !
+ pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
+ !
+ ! betaPE non-linear anomaly
+ zn2 = BPE002
+ !
+ zn1 = (BPE011)*zt &
+ & + BPE101*zs+BPE001
+ !
+ zn0 = (((BPE030)*zt &
+ & + BPE120*zs+BPE020)*zt &
+ & + (BPE210*zs+BPE110)*zs+BPE010)*zt &
+ & + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
+ !
+ zn = ( zn2 * zh + zn1 ) * zh + zn0
+ !
+ pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
+ !
+ END DO
+ END DO
+ END DO
+ !
+ CASE( np_seos ) !== Vallis (2006) simplified EOS ==!
+ !
+ DO jk = 1, jpkm1
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0)
+ zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
+ zh = gdept_n(ji,jj,jk) ! depth in meters at t-point
+ ztm = tmask(ji,jj,jk) ! tmask
+ zn = 0.5_wp * zh * r1_rau0 * ztm
+ ! ! Potential Energy
+ ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
+ ! ! alphaPE
+ pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
+ pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn
+ !
+ END DO
+ END DO
+ END DO
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,cform_err)
+ IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
+ nstop = nstop + 1
+ !
+ END SELECT
+ !
+ IF( ln_timing ) CALL timing_stop('eos_pen')
+ !
+ END SUBROUTINE eos_pen
+
+
+ SUBROUTINE eos_init
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE eos_init ***
+ !!
+ !! ** Purpose : initializations for the equation of state
+ !!
+ !! ** Method : Read the namelist nameos and control the parameters
+ !!----------------------------------------------------------------------
+ INTEGER :: ios ! local integer
+ INTEGER :: ioptio ! local integer
+ !!
+ NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1, &
+ & rn_lambda2, rn_mu2, rn_nu
+ !!----------------------------------------------------------------------
+ !
+ REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state
+ READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
+ !
+ REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state
+ READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
+902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
+ IF(lwm) WRITE( numond, nameos )
+ !
+ rau0 = 1026._wp !: volumic mass of reference [kg/m3]
+ rcp = 3991.86795711963_wp !: heat capacity [J/K]
+ !
+ IF(lwp) THEN ! Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'eos_init : equation of state'
+ WRITE(numout,*) '~~~~~~~~'
+ WRITE(numout,*) ' Namelist nameos : Chosen the Equation Of Seawater (EOS)'
+ WRITE(numout,*) ' TEOS-10 : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_TEOS10 = ', ln_TEOS10
+ WRITE(numout,*) ' EOS-80 : rho=F(Potential Temperature, Practical Salinity, depth) ln_EOS80 = ', ln_EOS80
+ WRITE(numout,*) ' S-EOS : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_SEOS = ', ln_SEOS
+ ENDIF
+
+ ! Check options for equation of state & set neos based on logical flags
+ ioptio = 0
+ IF( ln_TEOS10 ) THEN ; ioptio = ioptio+1 ; neos = np_teos10 ; ENDIF
+ IF( ln_EOS80 ) THEN ; ioptio = ioptio+1 ; neos = np_eos80 ; ENDIF
+ IF( ln_SEOS ) THEN ; ioptio = ioptio+1 ; neos = np_seos ; ENDIF
+ IF( ioptio /= 1 ) CALL ctl_stop("Exactly one equation of state option must be selected")
+ !
+ SELECT CASE( neos ) ! check option
+ !
+ CASE( np_teos10 ) !== polynomial TEOS-10 ==!
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' ==>>> use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
+ !
+ l_useCT = .TRUE. ! model temperature is Conservative temperature
+ !
+ rdeltaS = 32._wp
+ r1_S0 = 0.875_wp/35.16504_wp
+ r1_T0 = 1._wp/40._wp
+ r1_Z0 = 1.e-4_wp
+ !
+ EOS000 = 8.0189615746e+02_wp
+ EOS100 = 8.6672408165e+02_wp
+ EOS200 = -1.7864682637e+03_wp
+ EOS300 = 2.0375295546e+03_wp
+ EOS400 = -1.2849161071e+03_wp
+ EOS500 = 4.3227585684e+02_wp
+ EOS600 = -6.0579916612e+01_wp
+ EOS010 = 2.6010145068e+01_wp
+ EOS110 = -6.5281885265e+01_wp
+ EOS210 = 8.1770425108e+01_wp
+ EOS310 = -5.6888046321e+01_wp
+ EOS410 = 1.7681814114e+01_wp
+ EOS510 = -1.9193502195_wp
+ EOS020 = -3.7074170417e+01_wp
+ EOS120 = 6.1548258127e+01_wp
+ EOS220 = -6.0362551501e+01_wp
+ EOS320 = 2.9130021253e+01_wp
+ EOS420 = -5.4723692739_wp
+ EOS030 = 2.1661789529e+01_wp
+ EOS130 = -3.3449108469e+01_wp
+ EOS230 = 1.9717078466e+01_wp
+ EOS330 = -3.1742946532_wp
+ EOS040 = -8.3627885467_wp
+ EOS140 = 1.1311538584e+01_wp
+ EOS240 = -5.3563304045_wp
+ EOS050 = 5.4048723791e-01_wp
+ EOS150 = 4.8169980163e-01_wp
+ EOS060 = -1.9083568888e-01_wp
+ EOS001 = 1.9681925209e+01_wp
+ EOS101 = -4.2549998214e+01_wp
+ EOS201 = 5.0774768218e+01_wp
+ EOS301 = -3.0938076334e+01_wp
+ EOS401 = 6.6051753097_wp
+ EOS011 = -1.3336301113e+01_wp
+ EOS111 = -4.4870114575_wp
+ EOS211 = 5.0042598061_wp
+ EOS311 = -6.5399043664e-01_wp
+ EOS021 = 6.7080479603_wp
+ EOS121 = 3.5063081279_wp
+ EOS221 = -1.8795372996_wp
+ EOS031 = -2.4649669534_wp
+ EOS131 = -5.5077101279e-01_wp
+ EOS041 = 5.5927935970e-01_wp
+ EOS002 = 2.0660924175_wp
+ EOS102 = -4.9527603989_wp
+ EOS202 = 2.5019633244_wp
+ EOS012 = 2.0564311499_wp
+ EOS112 = -2.1311365518e-01_wp
+ EOS022 = -1.2419983026_wp
+ EOS003 = -2.3342758797e-02_wp
+ EOS103 = -1.8507636718e-02_wp
+ EOS013 = 3.7969820455e-01_wp
+ !
+ ALP000 = -6.5025362670e-01_wp
+ ALP100 = 1.6320471316_wp
+ ALP200 = -2.0442606277_wp
+ ALP300 = 1.4222011580_wp
+ ALP400 = -4.4204535284e-01_wp
+ ALP500 = 4.7983755487e-02_wp
+ ALP010 = 1.8537085209_wp
+ ALP110 = -3.0774129064_wp
+ ALP210 = 3.0181275751_wp
+ ALP310 = -1.4565010626_wp
+ ALP410 = 2.7361846370e-01_wp
+ ALP020 = -1.6246342147_wp
+ ALP120 = 2.5086831352_wp
+ ALP220 = -1.4787808849_wp
+ ALP320 = 2.3807209899e-01_wp
+ ALP030 = 8.3627885467e-01_wp
+ ALP130 = -1.1311538584_wp
+ ALP230 = 5.3563304045e-01_wp
+ ALP040 = -6.7560904739e-02_wp
+ ALP140 = -6.0212475204e-02_wp
+ ALP050 = 2.8625353333e-02_wp
+ ALP001 = 3.3340752782e-01_wp
+ ALP101 = 1.1217528644e-01_wp
+ ALP201 = -1.2510649515e-01_wp
+ ALP301 = 1.6349760916e-02_wp
+ ALP011 = -3.3540239802e-01_wp
+ ALP111 = -1.7531540640e-01_wp
+ ALP211 = 9.3976864981e-02_wp
+ ALP021 = 1.8487252150e-01_wp
+ ALP121 = 4.1307825959e-02_wp
+ ALP031 = -5.5927935970e-02_wp
+ ALP002 = -5.1410778748e-02_wp
+ ALP102 = 5.3278413794e-03_wp
+ ALP012 = 6.2099915132e-02_wp
+ ALP003 = -9.4924551138e-03_wp
+ !
+ BET000 = 1.0783203594e+01_wp
+ BET100 = -4.4452095908e+01_wp
+ BET200 = 7.6048755820e+01_wp
+ BET300 = -6.3944280668e+01_wp
+ BET400 = 2.6890441098e+01_wp
+ BET500 = -4.5221697773_wp
+ BET010 = -8.1219372432e-01_wp
+ BET110 = 2.0346663041_wp
+ BET210 = -2.1232895170_wp
+ BET310 = 8.7994140485e-01_wp
+ BET410 = -1.1939638360e-01_wp
+ BET020 = 7.6574242289e-01_wp
+ BET120 = -1.5019813020_wp
+ BET220 = 1.0872489522_wp
+ BET320 = -2.7233429080e-01_wp
+ BET030 = -4.1615152308e-01_wp
+ BET130 = 4.9061350869e-01_wp
+ BET230 = -1.1847737788e-01_wp
+ BET040 = 1.4073062708e-01_wp
+ BET140 = -1.3327978879e-01_wp
+ BET050 = 5.9929880134e-03_wp
+ BET001 = -5.2937873009e-01_wp
+ BET101 = 1.2634116779_wp
+ BET201 = -1.1547328025_wp
+ BET301 = 3.2870876279e-01_wp
+ BET011 = -5.5824407214e-02_wp
+ BET111 = 1.2451933313e-01_wp
+ BET211 = -2.4409539932e-02_wp
+ BET021 = 4.3623149752e-02_wp
+ BET121 = -4.6767901790e-02_wp
+ BET031 = -6.8523260060e-03_wp
+ BET002 = -6.1618945251e-02_wp
+ BET102 = 6.2255521644e-02_wp
+ BET012 = -2.6514181169e-03_wp
+ BET003 = -2.3025968587e-04_wp
+ !
+ PEN000 = -9.8409626043_wp
+ PEN100 = 2.1274999107e+01_wp
+ PEN200 = -2.5387384109e+01_wp
+ PEN300 = 1.5469038167e+01_wp
+ PEN400 = -3.3025876549_wp
+ PEN010 = 6.6681505563_wp
+ PEN110 = 2.2435057288_wp
+ PEN210 = -2.5021299030_wp
+ PEN310 = 3.2699521832e-01_wp
+ PEN020 = -3.3540239802_wp
+ PEN120 = -1.7531540640_wp
+ PEN220 = 9.3976864981e-01_wp
+ PEN030 = 1.2324834767_wp
+ PEN130 = 2.7538550639e-01_wp
+ PEN040 = -2.7963967985e-01_wp
+ PEN001 = -1.3773949450_wp
+ PEN101 = 3.3018402659_wp
+ PEN201 = -1.6679755496_wp
+ PEN011 = -1.3709540999_wp
+ PEN111 = 1.4207577012e-01_wp
+ PEN021 = 8.2799886843e-01_wp
+ PEN002 = 1.7507069098e-02_wp
+ PEN102 = 1.3880727538e-02_wp
+ PEN012 = -2.8477365341e-01_wp
+ !
+ APE000 = -1.6670376391e-01_wp
+ APE100 = -5.6087643219e-02_wp
+ APE200 = 6.2553247576e-02_wp
+ APE300 = -8.1748804580e-03_wp
+ APE010 = 1.6770119901e-01_wp
+ APE110 = 8.7657703198e-02_wp
+ APE210 = -4.6988432490e-02_wp
+ APE020 = -9.2436260751e-02_wp
+ APE120 = -2.0653912979e-02_wp
+ APE030 = 2.7963967985e-02_wp
+ APE001 = 3.4273852498e-02_wp
+ APE101 = -3.5518942529e-03_wp
+ APE011 = -4.1399943421e-02_wp
+ APE002 = 7.1193413354e-03_wp
+ !
+ BPE000 = 2.6468936504e-01_wp
+ BPE100 = -6.3170583896e-01_wp
+ BPE200 = 5.7736640125e-01_wp
+ BPE300 = -1.6435438140e-01_wp
+ BPE010 = 2.7912203607e-02_wp
+ BPE110 = -6.2259666565e-02_wp
+ BPE210 = 1.2204769966e-02_wp
+ BPE020 = -2.1811574876e-02_wp
+ BPE120 = 2.3383950895e-02_wp
+ BPE030 = 3.4261630030e-03_wp
+ BPE001 = 4.1079296834e-02_wp
+ BPE101 = -4.1503681096e-02_wp
+ BPE011 = 1.7676120780e-03_wp
+ BPE002 = 1.7269476440e-04_wp
+ !
+ CASE( np_eos80 ) !== polynomial EOS-80 formulation ==!
+ !
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' ==>>> use of EOS-80 equation of state (pot. temp. and pract. salinity)'
+ !
+ l_useCT = .FALSE. ! model temperature is Potential temperature
+ rdeltaS = 20._wp
+ r1_S0 = 1._wp/40._wp
+ r1_T0 = 1._wp/40._wp
+ r1_Z0 = 1.e-4_wp
+ !
+ EOS000 = 9.5356891948e+02_wp
+ EOS100 = 1.7136499189e+02_wp
+ EOS200 = -3.7501039454e+02_wp
+ EOS300 = 5.1856810420e+02_wp
+ EOS400 = -3.7264470465e+02_wp
+ EOS500 = 1.4302533998e+02_wp
+ EOS600 = -2.2856621162e+01_wp
+ EOS010 = 1.0087518651e+01_wp
+ EOS110 = -1.3647741861e+01_wp
+ EOS210 = 8.8478359933_wp
+ EOS310 = -7.2329388377_wp
+ EOS410 = 1.4774410611_wp
+ EOS510 = 2.0036720553e-01_wp
+ EOS020 = -2.5579830599e+01_wp
+ EOS120 = 2.4043512327e+01_wp
+ EOS220 = -1.6807503990e+01_wp
+ EOS320 = 8.3811577084_wp
+ EOS420 = -1.9771060192_wp
+ EOS030 = 1.6846451198e+01_wp
+ EOS130 = -2.1482926901e+01_wp
+ EOS230 = 1.0108954054e+01_wp
+ EOS330 = -6.2675951440e-01_wp
+ EOS040 = -8.0812310102_wp
+ EOS140 = 1.0102374985e+01_wp
+ EOS240 = -4.8340368631_wp
+ EOS050 = 1.2079167803_wp
+ EOS150 = 1.1515380987e-01_wp
+ EOS060 = -2.4520288837e-01_wp
+ EOS001 = 1.0748601068e+01_wp
+ EOS101 = -1.7817043500e+01_wp
+ EOS201 = 2.2181366768e+01_wp
+ EOS301 = -1.6750916338e+01_wp
+ EOS401 = 4.1202230403_wp
+ EOS011 = -1.5852644587e+01_wp
+ EOS111 = -7.6639383522e-01_wp
+ EOS211 = 4.1144627302_wp
+ EOS311 = -6.6955877448e-01_wp
+ EOS021 = 9.9994861860_wp
+ EOS121 = -1.9467067787e-01_wp
+ EOS221 = -1.2177554330_wp
+ EOS031 = -3.4866102017_wp
+ EOS131 = 2.2229155620e-01_wp
+ EOS041 = 5.9503008642e-01_wp
+ EOS002 = 1.0375676547_wp
+ EOS102 = -3.4249470629_wp
+ EOS202 = 2.0542026429_wp
+ EOS012 = 2.1836324814_wp
+ EOS112 = -3.4453674320e-01_wp
+ EOS022 = -1.2548163097_wp
+ EOS003 = 1.8729078427e-02_wp
+ EOS103 = -5.7238495240e-02_wp
+ EOS013 = 3.8306136687e-01_wp
+ !
+ ALP000 = -2.5218796628e-01_wp
+ ALP100 = 3.4119354654e-01_wp
+ ALP200 = -2.2119589983e-01_wp
+ ALP300 = 1.8082347094e-01_wp
+ ALP400 = -3.6936026529e-02_wp
+ ALP500 = -5.0091801383e-03_wp
+ ALP010 = 1.2789915300_wp
+ ALP110 = -1.2021756164_wp
+ ALP210 = 8.4037519952e-01_wp
+ ALP310 = -4.1905788542e-01_wp
+ ALP410 = 9.8855300959e-02_wp
+ ALP020 = -1.2634838399_wp
+ ALP120 = 1.6112195176_wp
+ ALP220 = -7.5817155402e-01_wp
+ ALP320 = 4.7006963580e-02_wp
+ ALP030 = 8.0812310102e-01_wp
+ ALP130 = -1.0102374985_wp
+ ALP230 = 4.8340368631e-01_wp
+ ALP040 = -1.5098959754e-01_wp
+ ALP140 = -1.4394226233e-02_wp
+ ALP050 = 3.6780433255e-02_wp
+ ALP001 = 3.9631611467e-01_wp
+ ALP101 = 1.9159845880e-02_wp
+ ALP201 = -1.0286156825e-01_wp
+ ALP301 = 1.6738969362e-02_wp
+ ALP011 = -4.9997430930e-01_wp
+ ALP111 = 9.7335338937e-03_wp
+ ALP211 = 6.0887771651e-02_wp
+ ALP021 = 2.6149576513e-01_wp
+ ALP121 = -1.6671866715e-02_wp
+ ALP031 = -5.9503008642e-02_wp
+ ALP002 = -5.4590812035e-02_wp
+ ALP102 = 8.6134185799e-03_wp
+ ALP012 = 6.2740815484e-02_wp
+ ALP003 = -9.5765341718e-03_wp
+ !
+ BET000 = 2.1420623987_wp
+ BET100 = -9.3752598635_wp
+ BET200 = 1.9446303907e+01_wp
+ BET300 = -1.8632235232e+01_wp
+ BET400 = 8.9390837485_wp
+ BET500 = -1.7142465871_wp
+ BET010 = -1.7059677327e-01_wp
+ BET110 = 2.2119589983e-01_wp
+ BET210 = -2.7123520642e-01_wp
+ BET310 = 7.3872053057e-02_wp
+ BET410 = 1.2522950346e-02_wp
+ BET020 = 3.0054390409e-01_wp
+ BET120 = -4.2018759976e-01_wp
+ BET220 = 3.1429341406e-01_wp
+ BET320 = -9.8855300959e-02_wp
+ BET030 = -2.6853658626e-01_wp
+ BET130 = 2.5272385134e-01_wp
+ BET230 = -2.3503481790e-02_wp
+ BET040 = 1.2627968731e-01_wp
+ BET140 = -1.2085092158e-01_wp
+ BET050 = 1.4394226233e-03_wp
+ BET001 = -2.2271304375e-01_wp
+ BET101 = 5.5453416919e-01_wp
+ BET201 = -6.2815936268e-01_wp
+ BET301 = 2.0601115202e-01_wp
+ BET011 = -9.5799229402e-03_wp
+ BET111 = 1.0286156825e-01_wp
+ BET211 = -2.5108454043e-02_wp
+ BET021 = -2.4333834734e-03_wp
+ BET121 = -3.0443885826e-02_wp
+ BET031 = 2.7786444526e-03_wp
+ BET002 = -4.2811838287e-02_wp
+ BET102 = 5.1355066072e-02_wp
+ BET012 = -4.3067092900e-03_wp
+ BET003 = -7.1548119050e-04_wp
+ !
+ PEN000 = -5.3743005340_wp
+ PEN100 = 8.9085217499_wp
+ PEN200 = -1.1090683384e+01_wp
+ PEN300 = 8.3754581690_wp
+ PEN400 = -2.0601115202_wp
+ PEN010 = 7.9263222935_wp
+ PEN110 = 3.8319691761e-01_wp
+ PEN210 = -2.0572313651_wp
+ PEN310 = 3.3477938724e-01_wp
+ PEN020 = -4.9997430930_wp
+ PEN120 = 9.7335338937e-02_wp
+ PEN220 = 6.0887771651e-01_wp
+ PEN030 = 1.7433051009_wp
+ PEN130 = -1.1114577810e-01_wp
+ PEN040 = -2.9751504321e-01_wp
+ PEN001 = -6.9171176978e-01_wp
+ PEN101 = 2.2832980419_wp
+ PEN201 = -1.3694684286_wp
+ PEN011 = -1.4557549876_wp
+ PEN111 = 2.2969116213e-01_wp
+ PEN021 = 8.3654420645e-01_wp
+ PEN002 = -1.4046808820e-02_wp
+ PEN102 = 4.2928871430e-02_wp
+ PEN012 = -2.8729602515e-01_wp
+ !
+ APE000 = -1.9815805734e-01_wp
+ APE100 = -9.5799229402e-03_wp
+ APE200 = 5.1430784127e-02_wp
+ APE300 = -8.3694846809e-03_wp
+ APE010 = 2.4998715465e-01_wp
+ APE110 = -4.8667669469e-03_wp
+ APE210 = -3.0443885826e-02_wp
+ APE020 = -1.3074788257e-01_wp
+ APE120 = 8.3359333577e-03_wp
+ APE030 = 2.9751504321e-02_wp
+ APE001 = 3.6393874690e-02_wp
+ APE101 = -5.7422790533e-03_wp
+ APE011 = -4.1827210323e-02_wp
+ APE002 = 7.1824006288e-03_wp
+ !
+ BPE000 = 1.1135652187e-01_wp
+ BPE100 = -2.7726708459e-01_wp
+ BPE200 = 3.1407968134e-01_wp
+ BPE300 = -1.0300557601e-01_wp
+ BPE010 = 4.7899614701e-03_wp
+ BPE110 = -5.1430784127e-02_wp
+ BPE210 = 1.2554227021e-02_wp
+ BPE020 = 1.2166917367e-03_wp
+ BPE120 = 1.5221942913e-02_wp
+ BPE030 = -1.3893222263e-03_wp
+ BPE001 = 2.8541225524e-02_wp
+ BPE101 = -3.4236710714e-02_wp
+ BPE011 = 2.8711395266e-03_wp
+ BPE002 = 5.3661089288e-04_wp
+ !
+ CASE( np_seos ) !== Simplified EOS ==!
+ IF(lwp) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) ' ==>>> use of simplified eos: '
+ WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
+ WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rau0'
+ WRITE(numout,*) ' with the following coefficients :'
+ WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0
+ WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0
+ WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1
+ WRITE(numout,*) ' cabbeling coef. rn_lambda2 = ', rn_lambda2
+ WRITE(numout,*) ' thermobar. coef. rn_mu1 = ', rn_mu1
+ WRITE(numout,*) ' thermobar. coef. rn_mu2 = ', rn_mu2
+ WRITE(numout,*) ' 2nd cabbel. coef. rn_nu = ', rn_nu
+ WRITE(numout,*) ' Caution: rn_beta0=0 incompatible with ddm parameterization '
+ ENDIF
+ l_useCT = .TRUE. ! Use conservative temperature
+ !
+ CASE DEFAULT !== ERROR in neos ==!
+ WRITE(ctmp1,*) ' bad flag value for neos = ', neos, '. You should never see this error'
+ CALL ctl_stop( ctmp1 )
+ !
+ END SELECT
+ !
+ rau0_rcp = rau0 * rcp
+ r1_rau0 = 1._wp / rau0
+ r1_rcp = 1._wp / rcp
+ r1_rau0_rcp = 1._wp / rau0_rcp
+ !
+ IF(lwp) THEN
+ IF( l_useCT ) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) ' ==>>> model uses Conservative Temperature'
+ WRITE(numout,*) ' Important: model must be initialized with CT and SA fields'
+ ELSE
+ WRITE(numout,*)
+ WRITE(numout,*) ' ==>>> model does not use Conservative Temperature'
+ ENDIF
+ ENDIF
+ !
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' Associated physical constant'
+ IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3'
+ IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg'
+ IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin'
+ IF(lwp) WRITE(numout,*) ' rau0 * rcp rau0_rcp = ', rau0_rcp
+ IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp
+ !
+ END SUBROUTINE eos_init
+
+ !!======================================================================
+END MODULE eosbn2
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/tide_mod.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/tide_mod.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/tide_mod.F90 (revision 10118)
@@ -0,0 +1,417 @@
+MODULE tide_mod
+ !!======================================================================
+ !! *** MODULE tide_mod ***
+ !! Compute nodal modulations corrections and pulsations
+ !!======================================================================
+ !! History : 1.0 ! 2007 (O. Le Galloudec) Original code
+ !!----------------------------------------------------------------------
+ USE dom_oce ! ocean space and time domain
+ USE phycst ! physical constant
+ USE daymod ! calendar
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC tide_harmo ! called by tideini and diaharm modules
+ PUBLIC tide_init_Wave ! called by tideini and diaharm modules
+
+ INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 19 !: maximum number of harmonic
+
+ TYPE, PUBLIC :: tide
+ CHARACTER(LEN=4) :: cname_tide
+ REAL(wp) :: equitide
+ INTEGER :: nutide
+ INTEGER :: nt, ns, nh, np, np1, shift
+ INTEGER :: nksi, nnu0, nnu1, nnu2, R
+ INTEGER :: nformula
+ END TYPE tide
+
+ TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave !:
+
+ REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1 ! astronomic angles
+ REAL(wp) :: sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R !
+ REAL(wp) :: sh_I, sh_x1ra, sh_N !
+
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id: tide_mod.F90 9598 2018-05-15 22:47:16Z nicolasmartin $
+ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE tide_init_Wave
+# include "tide.h90"
+ END SUBROUTINE tide_init_Wave
+
+
+ SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
+ !!----------------------------------------------------------------------
+ !!----------------------------------------------------------------------
+ INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
+ INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
+ REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
+ REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
+ !!----------------------------------------------------------------------
+ !
+ CALL astronomic_angle
+ CALL tide_pulse( pomega, ktide ,kc )
+
+! jc: Comodo IW case, remove nodal corrections:
+! CALL tide_vuf ( pvt, put, pcor, ktide ,kc )
+ pvt(:) = 0._wp
+ put(:) = 0._wp
+ pcor(:) = 1._wp
+! jc
+
+ END SUBROUTINE tide_harmo
+
+
+ SUBROUTINE astronomic_angle
+ !!----------------------------------------------------------------------
+ !! tj is time elapsed since 1st January 1900, 0 hour, counted in julian
+ !! century (e.g. time in days divide by 36525)
+ !!----------------------------------------------------------------------
+ REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
+ REAL(wp) :: zqy , zsy, zday, zdj, zhfrac
+ !!----------------------------------------------------------------------
+ !
+ zqy = AINT( (nyear-1901.)/4. )
+ zsy = nyear - 1900.
+ !
+ zdj = dayjul( nyear, nmonth, nday )
+ zday = zdj + zqy - 1.
+ !
+ zhfrac = nsec_day / 3600.
+ !
+ !----------------------------------------------------------------------
+ ! Sh_n Longitude of ascending lunar node
+ !----------------------------------------------------------------------
+ sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
+ !----------------------------------------------------------------------
+ ! T mean solar angle (Greenwhich time)
+ !----------------------------------------------------------------------
+ sh_T=(180.+zhfrac*(360./24.))*rad
+ !----------------------------------------------------------------------
+ ! h mean solar Longitude
+ !----------------------------------------------------------------------
+ sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
+ !----------------------------------------------------------------------
+ ! s mean lunar Longitude
+ !----------------------------------------------------------------------
+ sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
+ !----------------------------------------------------------------------
+ ! p1 Longitude of solar perigee
+ !----------------------------------------------------------------------
+ sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
+ !----------------------------------------------------------------------
+ ! p Longitude of lunar perigee
+ !----------------------------------------------------------------------
+ sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
+
+ sh_N = MOD( sh_N ,2*rpi )
+ sh_s = MOD( sh_s ,2*rpi )
+ sh_h = MOD( sh_h, 2*rpi )
+ sh_p = MOD( sh_p, 2*rpi )
+ sh_p1= MOD( sh_p1,2*rpi )
+
+ cosI = 0.913694997 -0.035692561 *cos(sh_N)
+
+ sh_I = ACOS( cosI )
+
+ sin2I = sin(sh_I)
+ sh_tgn2 = tan(sh_N/2.0)
+
+ at1=atan(1.01883*sh_tgn2)
+ at2=atan(0.64412*sh_tgn2)
+
+ sh_xi=-at1-at2+sh_N
+
+ IF( sh_N > rpi ) sh_xi=sh_xi-2.0*rpi
+
+ sh_nu = at1 - at2
+
+ !----------------------------------------------------------------------
+ ! For constituents l2 k1 k2
+ !----------------------------------------------------------------------
+
+ tgI2 = tan(sh_I/2.0)
+ P1 = sh_p-sh_xi
+
+ t2 = tgI2*tgI2
+ t4 = t2*t2
+ sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
+
+ p = sin(2.0*P1)
+ q = 1.0/(6.0*t2)-cos(2.0*P1)
+ sh_R = atan(p/q)
+
+ p = sin(2.0*sh_I)*sin(sh_nu)
+ q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
+ sh_nuprim = atan(p/q)
+
+ s2 = sin(sh_I)*sin(sh_I)
+ p = s2*sin(2.0*sh_nu)
+ q = s2*cos(2.0*sh_nu)+0.0727
+ sh_nusec = 0.5*atan(p/q)
+ !
+ END SUBROUTINE astronomic_angle
+
+
+ SUBROUTINE tide_pulse( pomega, ktide ,kc )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tide_pulse ***
+ !!
+ !! ** Purpose : Compute tidal frequencies
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
+ INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
+ REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
+ !
+ INTEGER :: jh
+ REAL(wp) :: zscale
+ REAL(wp) :: zomega_T = 13149000.0_wp
+ REAL(wp) :: zomega_s = 481267.892_wp
+ REAL(wp) :: zomega_h = 36000.76892_wp
+ REAL(wp) :: zomega_p = 4069.0322056_wp
+ REAL(wp) :: zomega_n = 1934.1423972_wp
+ REAL(wp) :: zomega_p1= 1.719175_wp
+ !!----------------------------------------------------------------------
+ !
+ zscale = rad / ( 36525._wp * 86400._wp )
+ !
+ DO jh = 1, kc
+ pomega(jh) = ( zomega_T * Wave( ktide(jh) )%nT &
+ & + zomega_s * Wave( ktide(jh) )%ns &
+ & + zomega_h * Wave( ktide(jh) )%nh &
+ & + zomega_p * Wave( ktide(jh) )%np &
+ & + zomega_p1* Wave( ktide(jh) )%np1 ) * zscale
+ END DO
+ !
+ END SUBROUTINE tide_pulse
+
+
+ SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tide_vuf ***
+ !!
+ !! ** Purpose : Compute nodal modulation corrections
+ !!
+ !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
+ !! ut: Phase correction u due to nodal motion (radians)
+ !! ft: Nodal correction factor
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
+ INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
+ REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
+ !
+ INTEGER :: jh ! dummy loop index
+ !!----------------------------------------------------------------------
+ !
+ DO jh = 1, kc
+ ! Phase of the tidal potential relative to the Greenwhich
+ ! meridian (e.g. the position of the fictuous celestial body). Units are radian:
+ pvt(jh) = sh_T * Wave( ktide(jh) )%nT &
+ & + sh_s * Wave( ktide(jh) )%ns &
+ & + sh_h * Wave( ktide(jh) )%nh &
+ & + sh_p * Wave( ktide(jh) )%np &
+ & + sh_p1* Wave( ktide(jh) )%np1 &
+ & + Wave( ktide(jh) )%shift * rad
+ !
+ ! Phase correction u due to nodal motion. Units are radian:
+ put(jh) = sh_xi * Wave( ktide(jh) )%nksi &
+ & + sh_nu * Wave( ktide(jh) )%nnu0 &
+ & + sh_nuprim * Wave( ktide(jh) )%nnu1 &
+ & + sh_nusec * Wave( ktide(jh) )%nnu2 &
+ & + sh_R * Wave( ktide(jh) )%R
+
+ ! Nodal correction factor:
+ pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
+ END DO
+ !
+ END SUBROUTINE tide_vuf
+
+
+ RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
+ !!----------------------------------------------------------------------
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT(in) :: kformula
+ !
+ REAL(wp) :: zf
+ REAL(wp) :: zs, zf1, zf2
+ !!----------------------------------------------------------------------
+ !
+ SELECT CASE( kformula )
+ !
+ CASE( 0 ) !== formule 0, solar waves
+ zf = 1.0
+ !
+ CASE( 1 ) !== formule 1, compound waves (78 x 78)
+ zf=nodal_factort(78)
+ zf = zf * zf
+ !
+ CASE ( 2 ) !== formule 2, compound waves (78 x 0) === (78)
+ zf1= nodal_factort(78)
+ zf = nodal_factort( 0)
+ zf = zf1 * zf
+ !
+ CASE ( 4 ) !== formule 4, compound waves (78 x 235)
+ zf1 = nodal_factort( 78)
+ zf = nodal_factort(235)
+ zf = zf1 * zf
+ !
+ CASE ( 5 ) !== formule 5, compound waves (78 *78 x 235)
+ zf1 = nodal_factort( 78)
+ zf = nodal_factort(235)
+ zf = zf * zf1 * zf1
+ !
+ CASE ( 6 ) !== formule 6, compound waves (78 *78 x 0)
+ zf1 = nodal_factort(78)
+ zf = nodal_factort( 0)
+ zf = zf * zf1 * zf1
+ !
+ CASE( 7 ) !== formule 7, compound waves (75 x 75)
+ zf = nodal_factort(75)
+ zf = zf * zf
+ !
+ CASE( 8 ) !== formule 8, compound waves (78 x 0 x 235)
+ zf = nodal_factort( 78)
+ zf1 = nodal_factort( 0)
+ zf2 = nodal_factort(235)
+ zf = zf * zf1 * zf2
+ !
+ CASE( 9 ) !== formule 9, compound waves (78 x 0 x 227)
+ zf = nodal_factort( 78)
+ zf1 = nodal_factort( 0)
+ zf2 = nodal_factort(227)
+ zf = zf * zf1 * zf2
+ !
+ CASE( 10 ) !== formule 10, compound waves (78 x 227)
+ zf = nodal_factort( 78)
+ zf1 = nodal_factort(227)
+ zf = zf * zf1
+ !
+ CASE( 11 ) !== formule 11, compound waves (75 x 0)
+!!gm bug???? zf 2 fois !
+ zf = nodal_factort(75)
+ zf1 = nodal_factort( 0)
+ zf = zf * zf1
+ !
+ CASE( 12 ) !== formule 12, compound waves (78 x 78 x 78 x 0)
+ zf1 = nodal_factort(78)
+ zf = nodal_factort( 0)
+ zf = zf * zf1 * zf1 * zf1
+ !
+ CASE( 13 ) !== formule 13, compound waves (78 x 75)
+ zf1 = nodal_factort(78)
+ zf = nodal_factort(75)
+ zf = zf * zf1
+ !
+ CASE( 14 ) !== formule 14, compound waves (235 x 0) === (235)
+ zf = nodal_factort(235)
+ zf1 = nodal_factort( 0)
+ zf = zf * zf1
+ !
+ CASE( 15 ) !== formule 15, compound waves (235 x 75)
+ zf = nodal_factort(235)
+ zf1 = nodal_factort( 75)
+ zf = zf * zf1
+ !
+ CASE( 16 ) !== formule 16, compound waves (78 x 0 x 0) === (78)
+ zf = nodal_factort(78)
+ zf1 = nodal_factort( 0)
+ zf = zf * zf1 * zf1
+ !
+ CASE( 17 ) !== formule 17, compound waves (227 x 0)
+ zf1 = nodal_factort(227)
+ zf = nodal_factort( 0)
+ zf = zf * zf1
+ !
+ CASE( 18 ) !== formule 18, compound waves (78 x 78 x 78 )
+ zf1 = nodal_factort(78)
+ zf = zf1 * zf1 * zf1
+ !
+ CASE( 19 ) !== formule 19, compound waves (78 x 0 x 0 x 0) === (78)
+!!gm bug2 ==>>> here identical to formule 16, a third multiplication by zf1 is missing
+ zf = nodal_factort(78)
+ zf1 = nodal_factort( 0)
+ zf = zf * zf1 * zf1
+ !
+ CASE( 73 ) !== formule 73
+ zs = sin(sh_I)
+ zf = (2./3.-zs*zs)/0.5021
+ !
+ CASE( 74 ) !== formule 74
+ zs = sin(sh_I)
+ zf = zs * zs / 0.1578
+ !
+ CASE( 75 ) !== formule 75
+ zs = cos(sh_I/2)
+ zf = sin(sh_I) * zs * zs / 0.3800
+ !
+ CASE( 76 ) !== formule 76
+ zf = sin(2*sh_I) / 0.7214
+ !
+ CASE( 77 ) !== formule 77
+ zs = sin(sh_I/2)
+ zf = sin(sh_I) * zs * zs / 0.0164
+ !
+ CASE( 78 ) !== formule 78
+ zs = cos(sh_I/2)
+ zf = zs * zs * zs * zs / 0.9154
+ !
+ CASE( 79 ) !== formule 79
+ zs = sin(sh_I)
+ zf = zs * zs / 0.1565
+ !
+ CASE( 144 ) !== formule 144
+ zs = sin(sh_I/2)
+ zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
+ !
+ CASE( 149 ) !== formule 149
+ zs = cos(sh_I/2)
+ zf = zs*zs*zs*zs*zs*zs / 0.8758
+ !
+ CASE( 215 ) !== formule 215
+ zs = cos(sh_I/2)
+ zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
+ !
+ CASE( 227 ) !== formule 227
+ zs = sin(2*sh_I)
+ zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
+ !
+ CASE ( 235 ) !== formule 235
+ zs = sin(sh_I)
+ zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
+ !
+ END SELECT
+ !
+ END FUNCTION nodal_factort
+
+
+ FUNCTION dayjul( kyr, kmonth, kday )
+ !!----------------------------------------------------------------------
+ !! *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
+ !!----------------------------------------------------------------------
+ INTEGER,INTENT(in) :: kyr, kmonth, kday
+ !
+ INTEGER,DIMENSION(12) :: idayt, idays
+ INTEGER :: inc, ji
+ REAL(wp) :: dayjul, zyq
+ !
+ DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
+ !!----------------------------------------------------------------------
+ !
+ idays(1) = 0.
+ idays(2) = 31.
+ inc = 0.
+ zyq = MOD( kyr-1900. , 4. )
+ IF( zyq == 0.) inc = 1.
+ DO ji = 3, 12
+ idays(ji)=idayt(ji)+inc
+ END DO
+ dayjul = idays(kmonth) + kday
+ !
+ END FUNCTION dayjul
+
+ !!======================================================================
+END MODULE tide_mod
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_fmask.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_fmask.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_fmask.F90 (revision 10118)
@@ -0,0 +1,150 @@
+MODULE usrdef_fmask
+ !!======================================================================
+ !! *** MODULE usrdef_fmask ***
+ !!
+ !! === ORCA configuration ===
+ !! (2 and 1 degrees)
+ !!
+ !! User defined : alteration of land/sea f-point mask in some straits
+ !!======================================================================
+ !! History : 4.0 ! 2016-06 (G. Madec, S. Flavoni) Original code
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! usr_def_fmask : alteration of f-point land/ocean mask in some straits
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers
+ USE dom_oce ! ocean space and time domain
+ !
+ USE in_out_manager ! I/O manager
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_mpp ! Massively Parallel Processing library
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC usr_def_fmask ! routine called by dommsk.F90
+
+ !! * Substitutions
+# include "vectopt_loop_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id:$
+ !! Software governed by the CeCILL licence (./LICENSE)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE usr_def_fmask( cd_cfg, kcfg, pfmsk )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE dom_msk ***
+ !!
+ !! ** Purpose : User defined alteration of the lateral boundary
+ !! condition on velocity.
+ !!
+ !! ** Method : Local change of the value of fmask at lateral ocean/land
+ !! boundary in straits in order to increase the viscous
+ !! boundary layer and thus reduce the transport through the
+ !! corresponding straits.
+ !! Here only alterations in ORCA R2 and R1 cases
+ !!
+ !! ** Action : fmask : land/ocean mask at f-point with increased value
+ !! in some user defined straits
+ !!----------------------------------------------------------------------
+ CHARACTER(len=*) , INTENT(in ) :: cd_cfg ! configuration name
+ INTEGER , INTENT(in ) :: kcfg ! configuration identifier
+ REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfmsk ! Ocean/Land f-point mask including lateral boundary cond.
+ !
+ INTEGER :: iif, iil, ii0, ii1, ii ! local integers
+ INTEGER :: ijf, ijl, ij0, ij1 ! - -
+ INTEGER :: isrow ! index for ORCA1 starting row
+ !!----------------------------------------------------------------------
+ !
+ IF( TRIM( cd_cfg ) == "orca" ) THEN !== ORCA Configurations ==!
+ !
+ SELECT CASE ( kcfg )
+ !
+ CASE( 2 ) ! R2 case
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R2: increase lateral friction near the following straits:'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+ !
+ IF(lwp) WRITE(numout,*) ' Gibraltar '
+ ij0 = 101 ; ij1 = 101 ! Gibraltar strait : partial slip (pfmsk=0.5)
+ ii0 = 139 ; ii1 = 140 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
+ ij0 = 102 ; ij1 = 102
+ ii0 = 139 ; ii1 = 140 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
+ !
+ IF(lwp) WRITE(numout,*) ' Bab el Mandeb '
+ ij0 = 87 ; ij1 = 88 ! Bab el Mandeb : partial slip (pfmsk=1)
+ ii0 = 160 ; ii1 = 160 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
+ ij0 = 88 ; ij1 = 88
+ ii0 = 159 ; ii1 = 159 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
+ !
+ ! We keep this as an example but it is instable in this case
+ !IF(lwp) WRITE(numout,*) ' Danish straits '
+ ! ij0 = 115 ; ij1 = 115 ! Danish straits : strong slip (pfmsk > 2)
+ ! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
+ ! ij0 = 116 ; ij1 = 116
+ ! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
+ !
+ CASE( 1 ) ! R1 case
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R1: increase lateral friction near the following straits:'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+!!gm ! This dirty section will be suppressed by simplification process:
+!!gm ! all this will come back in input files
+!!gm ! Currently these hard-wired indices relate to configuration with extend grid (jpjglo=332)
+ !
+ isrow = 332 - jpjglo
+ !
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) ' orca_r1: increase friction near the following straits : '
+ IF(lwp) WRITE(numout,*) ' Gibraltar '
+ ii0 = 282 ; ii1 = 283 ! Gibraltar Strait
+ ij0 = 241 - isrow ; ij1 = 241 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
+ !
+ IF(lwp) WRITE(numout,*) ' Bhosporus '
+ ii0 = 314 ; ii1 = 315 ! Bhosporus Strait
+ ij0 = 248 - isrow ; ij1 = 248 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
+ !
+ IF(lwp) WRITE(numout,*) ' Makassar (Top) '
+ ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top)
+ ij0 = 189 - isrow ; ij1 = 190 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
+ !
+ IF(lwp) WRITE(numout,*) ' Lombok '
+ ii0 = 44 ; ii1 = 44 ! Lombok Strait
+ ij0 = 164 - isrow ; ij1 = 165 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
+ !
+ IF(lwp) WRITE(numout,*) ' Ombai '
+ ii0 = 53 ; ii1 = 53 ! Ombai Strait
+ ij0 = 164 - isrow ; ij1 = 165 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
+ !
+ IF(lwp) WRITE(numout,*) ' Timor Passage '
+ ii0 = 56 ; ii1 = 56 ! Timor Passage
+ ij0 = 164 - isrow ; ij1 = 165 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
+ !
+ IF(lwp) WRITE(numout,*) ' West Halmahera '
+ ii0 = 58 ; ii1 = 58 ! West Halmahera Strait
+ ij0 = 181 - isrow ; ij1 = 182 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
+ !
+ IF(lwp) WRITE(numout,*) ' East Halmahera '
+ ii0 = 55 ; ii1 = 55 ! East Halmahera Strait
+ ij0 = 181 - isrow ; ij1 = 182 - isrow ; pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
+ !
+ CASE DEFAULT
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R', kcfg,' : NO alteration of fmask in specific straits '
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+ END SELECT
+ ELSE
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'usr_def_fmask : NO alteration of fmask in specific straits '
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+ ENDIF
+ !
+ CALL lbc_lnk( pfmsk, 'F', 1._wp ) ! Lateral boundary conditions on fmask
+ !
+ END SUBROUTINE usr_def_fmask
+
+ !!======================================================================
+END MODULE usrdef_fmask
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_istate.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_istate.F90 (revision 10118)
@@ -0,0 +1,124 @@
+MODULE usrdef_istate
+ !!======================================================================
+ !! *** MODULE usrdef_istate ***
+ !!
+ !! === COMODO-IW configuration ===
+ !!
+ !! User defined : set the initial state of a user configuration
+ !!======================================================================
+ !! History : 4.0 ! 2018-09 (J. Chanut) Original code
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! usr_def_istate : initial state in Temperature and salinity
+ !!----------------------------------------------------------------------
+ USE par_oce ! ocean space and time domain
+ USE phycst ! physical constants
+ USE eosbn2, ONLY : rn_a0
+ !
+ USE in_out_manager ! I/O manager
+ USE lib_mpp ! MPP library
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC usr_def_istate ! called in istate.F90
+
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id:$
+ !! Software governed by the CeCILL licence (./LICENSE)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE usr_def_istate ***
+ !!
+ !! ** Purpose : Initialization of the dynamics and tracers
+ !! for "COMODO internal waves" setup
+ !!
+ !! ** Method : - set temprature field
+ !! - set salinity field
+ !!----------------------------------------------------------------------
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pdept ! depth of t-point [m]
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: ptmask ! t-point ocean mask [m]
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pts ! T & S fields [Celsius ; g/kg]
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pu ! i-component of the velocity [m/s]
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pv ! j-component of the velocity [m/s]
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pssh ! sea-surface height
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ INTEGER :: kcase
+ REAL(wp):: znn ! Stratification
+ !!----------------------------------------------------------------------
+ !
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'usr_def_istate : analytical definition of initial state '
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ Ocean at rest, with an horizontally uniform T and S profiles'
+ !
+ pu (:,:,:) = 0._wp ! ocean at rest
+ pv (:,:,:) = 0._wp
+ pssh(:,:) = 0._wp
+ !
+ znn = 0.002_wp ! N [s-1]
+ kcase = 1
+ !
+ SELECT CASE ( kcase )
+ !
+ CASE( 1 )
+ DO jk = 1, jpk ! Linearly stratified
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1026.2_wp-rau0) / rn_a0 &
+ & -(znn**2) * rau0 /grav / rn_a0 &
+ & * (pdept(ji,jj,jk)-50._wp) &
+ & ) * ptmask(ji,jj,jk)
+ pts(ji,jj,jk,jp_sal) = 35._wp * ptmask(ji,jj,jk)
+ END DO
+ END DO
+ END DO
+ !
+ CASE( 2 )
+ DO jk=1, jpkm1 ! Interfacial waves:
+ DO jj=1,jpj ! Set a two layer system with
+ DO ji=1,jpi ! a density jump = 1.2 kg/m3 at z=50m
+ IF ( 0.5_wp*(pdept(ji,jj,jk)+pdept(ji,jj,jk+1))<= 50._wp ) THEN
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1025.8_wp-rau0) / rn_a0 ) * ptmask(ji,jj,jk)
+ ELSE
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1027.0_wp-rau0) / rn_a0 ) * ptmask(ji,jj,jk)
+ ENDIF
+ pts(ji,jj,jk,jp_sal) = 35._wp * ptmask(ji,jj,jk)
+ END DO
+ END DO
+ END DO
+ !
+ CASE( 3 )
+ DO jk=1, jpkm1 ! Interfacial waves:
+ DO jj=1,jpj ! stratification + mixed layer and thermocline
+ DO ji=1,jpi
+ IF (pdept(ji,jj,jk)< 50._wp) THEN
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1025.0_wp-rau0) / rn_a0 ) * ptmask(ji,jj,jk)
+ ELSEIF ((pdept(ji,jj,jk)>= 50._wp).AND.(pdept(ji,jj,jk)< 60._wp)) THEN
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1026.2_wp-rau0) / rn_a0 &
+ & -(znn**2) * rau0 /grav / rn_a0 &
+ & * (pdept(ji,jj,jk)-50._wp) &
+ & -1.2_wp * (pdept(ji,jj,jk)-60._wp) / 10._wp / rn_a0 &
+ & ) * ptmask(ji,jj,jk)
+ ELSEIF (pdept(ji,jj,jk)>= 60._wp) THEN
+ pts(ji,jj,jk,jp_tem) = (10._wp - (1026.2_wp-rau0) / rn_a0 &
+ & -(znn**2) * rau0 /grav / rn_a0 &
+ & * (pdept(ji,jj,jk)-50._wp) &
+ & ) * ptmask(ji,jj,jk)
+ ENDIF
+ pts(ji,jj,jk,jp_sal) = 35._wp * ptmask(ji,jj,jk)
+ END DO
+ END DO
+ END DO
+ !
+ END SELECT
+ !
+ END SUBROUTINE usr_def_istate
+
+ !!======================================================================
+END MODULE usrdef_istate
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_sbc.F90 (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/usrdef_sbc.F90 (revision 10118)
@@ -0,0 +1,86 @@
+MODULE usrdef_sbc
+ !!======================================================================
+ !! *** MODULE usrdef_sbc ***
+ !!
+ !! === COMODO-IW configuration ===
+ !!
+ !! User defined : surface forcing of a user configuration
+ !!======================================================================
+ !! History : 4.0 ! 2018-09 (J.Chanut) user defined interface
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! usr_def_sbc : user defined surface bounday conditions in OVERFLOW case
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers
+ USE dom_oce ! ocean space and time domain
+ USE sbc_oce ! Surface boundary condition: ocean fields
+ USE phycst ! physical constants
+ !
+ USE in_out_manager ! I/O manager
+ USE lib_mpp ! distribued memory computing library
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC usrdef_sbc_oce ! routine called in sbcmod module
+ PUBLIC usrdef_sbc_ice_tau ! routine called by icestp.F90 for ice dynamics
+ PUBLIC usrdef_sbc_ice_flx ! routine called by icestp.F90 for ice thermo
+
+ !! * Substitutions
+# include "vectopt_loop_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+ !! $Id$
+ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE usrdef_sbc_oce( kt )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE usr_def_sbc ***
+ !!
+ !! ** Purpose : provide at each time-step the surface boundary
+ !! condition, i.e. the momentum, heat and freshwater fluxes.
+ !!
+ !! ** Method : all 0 fields, for COMODO-IW case
+ !! CAUTION : never mask the surface stress field !
+ !!
+ !! ** Action : - set to ZERO all the ocean surface boundary condition, i.e.
+ !! utau, vtau, taum, wndm, qns, qsr, emp, sfx
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT(in) :: kt ! ocean time step
+ !!---------------------------------------------------------------------
+ !
+ IF( kt == nit000 ) THEN
+ !
+ IF(lwp) WRITE(numout,*)' usr_sbc : COMODO-IW case: NO surface forcing'
+ IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~ utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'
+ !
+ utau(:,:) = 0._wp
+ vtau(:,:) = 0._wp
+ taum(:,:) = 0._wp
+ wndm(:,:) = 0._wp
+ !
+ emp (:,:) = 0._wp
+ sfx (:,:) = 0._wp
+ qns (:,:) = 0._wp
+ qsr (:,:) = 0._wp
+ !
+ ENDIF
+ !
+ END SUBROUTINE usrdef_sbc_oce
+
+ SUBROUTINE usrdef_sbc_ice_tau( kt )
+ INTEGER, INTENT(in) :: kt ! ocean time step
+ END SUBROUTINE usrdef_sbc_ice_tau
+
+ SUBROUTINE usrdef_sbc_ice_flx( kt )
+ INTEGER, INTENT(in) :: kt ! ocean time step
+ END SUBROUTINE usrdef_sbc_ice_flx
+
+ !!======================================================================
+END MODULE usrdef_sbc
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/cpp_COMODO-IW.fcm
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/cpp_COMODO-IW.fcm (revision 10118)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/cpp_COMODO-IW.fcm (revision 10118)
@@ -0,0 +1,1 @@
+ bld::tool::fppkeys key_iomput key_mpp_mpi key_nosignedzero
Index: /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/demo_cfgs.txt
===================================================================
--- /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/demo_cfgs.txt (revision 10117)
+++ /NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/demo_cfgs.txt (revision 10118)
@@ -6,2 +6,3 @@
VORTEX OCE NST
WAD OCE
+COMODO-IW OCE