New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

ldftra.F90 in NEMO/trunk/src/OCE/LDF – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/LDF/ldftra.F90 @ 12276

Last change on this file since 12276 was 12276, checked in by cetlod, 4 years ago
trunk : merge in some cmip6 diagnostics into the trunk before copying it to release-4.0.2(-head). SETTE tests are OK and the is no difference with the revision 12248
Property svn:keywords set to `Id`
File size: 53.1 KB

Rev	Line
[3]	1	MODULE ldftra
	2	!!======================================================================
	3	!! * MODULE ldftra *
[5836]	4	!! Ocean physics: lateral diffusivity coefficients
[3]	5	!!=====================================================================
[5836]	6	!! History : ! 1997-07 (G. Madec) from inimix.F split in 2 routines
	7	!! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module
	8	!! 2.0 ! 2005-11 (G. Madec)
	9	!! 3.7 ! 2013-12 (F. Lemarie, G. Madec) restructuration/simplification of aht/aeiv specification,
	10	!! ! add velocity dependent coefficient and optional read in file
[3]	11	!!----------------------------------------------------------------------
[1601]	12
	13	!!----------------------------------------------------------------------
[3]	14	!! ldf_tra_init : initialization, namelist read, and parameters control
[5836]	15	!! ldf_tra : update lateral eddy diffusivity coefficients at each time step
	16	!! ldf_eiv_init : initialization of the eiv coeff. from namelist choices
	17	!! ldf_eiv : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability)
	18	!! ldf_eiv_trp : add to the input ocean transport the contribution of the EIV parametrization
	19	!! ldf_eiv_dia : diagnose the eddy induced velocity from the eiv streamfunction
[3]	20	!!----------------------------------------------------------------------
	21	USE oce ! ocean dynamics and tracers
	22	USE dom_oce ! ocean space and time domain
	23	USE phycst ! physical constants
[5836]	24	USE ldfslp ! lateral diffusion: slope of iso-neutral surfaces
	25	USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases
[7646]	26	USE diaptr
[5836]	27	!
[3]	28	USE in_out_manager ! I/O manager
[5836]	29	USE iom ! I/O module for ehanced bottom friction file
[3]	30	USE lib_mpp ! distribued memory computing library
	31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
[5836]	36	PUBLIC ldf_tra_init ! called by nemogcm.F90
	37	PUBLIC ldf_tra ! called by step.F90
	38	PUBLIC ldf_eiv_init ! called by nemogcm.F90
	39	PUBLIC ldf_eiv ! called by step.F90
	40	PUBLIC ldf_eiv_trp ! called by traadv.F90
	41	PUBLIC ldf_eiv_dia ! called by traldf_iso and traldf_iso_triad.F90
	42
	43	! !!* Namelist namtra_ldf : lateral mixing on tracers *
	44	! != Operator type =!
[9526]	45	LOGICAL , PUBLIC :: ln_traldf_OFF !: no operator: No explicit diffusion
[5836]	46	LOGICAL , PUBLIC :: ln_traldf_lap !: laplacian operator
	47	LOGICAL , PUBLIC :: ln_traldf_blp !: bilaplacian operator
	48	! != Direction of action =!
	49	LOGICAL , PUBLIC :: ln_traldf_lev !: iso-level direction
	50	LOGICAL , PUBLIC :: ln_traldf_hor !: horizontal (geopotential) direction
	51	! LOGICAL , PUBLIC :: ln_traldf_iso !: iso-neutral direction (see ldfslp)
[9490]	52	! != iso-neutral options =!
[5836]	53	! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp)
	54	LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction
	55	! LOGICAL , PUBLIC :: ln_triad_iso !: pure horizontal mixing in ML (see ldfslp)
	56	! LOGICAL , PUBLIC :: ln_botmix_triad !: mixing on bottom (see ldfslp)
	57	! REAL(wp), PUBLIC :: rn_sw_triad !: =1/0 switching triad / all 4 triads used (see ldfslp)
	58	! REAL(wp), PUBLIC :: rn_slpmax !: slope limit (see ldfslp)
	59	! != Coefficients =!
	60	INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef.
[9490]	61	! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case)
	62	! ! bht_0 = 1/12 Ud*Ld^3 (blp case)
	63	REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s]
	64	REAL(wp), PUBLIC :: rn_Ld !: lateral diffusive length [m]
[3]	65
[9490]	66	! !!* Namelist namtra_eiv : eddy induced velocity param. *
[5836]	67	! != Use/diagnose eiv =!
	68	LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag
	69	LOGICAL , PUBLIC :: ln_ldfeiv_dia !: diagnose & output eiv streamfunction and velocity (IOM)
	70	! != Coefficients =!
	71	INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff.
[9490]	72	REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s]
	73	REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m]
[5836]	74
[9490]	75	! ! Flag to control the type of lateral diffusive operator
	76	INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion
	77	INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend)
	78	! !! laplacian ! bilaplacian !
	79	INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator
	80	INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator
	81	INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator
	82
	83	INTEGER , PUBLIC :: nldf_tra = 0 !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals)
[5836]	84	LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef.
[9490]	85	LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef.
[5836]	86
	87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s]
	88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s]
	89
[9490]	90	REAL(wp) :: aht0, aei0 ! constant eddy coefficients (deduced from namelist values) [m2/s]
	91	REAL(wp) :: r1_2 = 0.5_wp ! =1/2
[5836]	92	REAL(wp) :: r1_4 = 0.25_wp ! =1/4
	93	REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12
	94
[3]	95	!! * Substitutions
	96	# include "vectopt_loop_substitute.h90"
[1601]	97	!!----------------------------------------------------------------------
[9598]	98	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	99	!! $Id$
[10068]	100	!! Software governed by the CeCILL license (see ./LICENSE)
[1601]	101	!!----------------------------------------------------------------------
[3]	102	CONTAINS
	103
	104	SUBROUTINE ldf_tra_init
	105	!!----------------------------------------------------------------------
	106	!! * ROUTINE ldf_tra_init *
	107	!!
[461]	108	!! ** Purpose : initializations of the tracer lateral mixing coeff.
[3]	109	!!
[5836]	110	!! ** Method : * the eddy diffusivity coef. specification depends on:
[3]	111	!!
[5836]	112	!! ln_traldf_lap = T laplacian operator
	113	!! ln_traldf_blp = T bilaplacian operator
	114	!!
	115	!! nn_aht_ijk_t = 0 => = constant
	116	!! !
	117	!! = 10 => = F(z) : constant with a reduction of 1/4 with depth
	118	!! !
	119	!! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file
	120	!! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
	121	!! = 21 = F(i,j,t) = F(growth rate of baroclinic instability)
	122	!! !
	123	!! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file
	124	!! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10)
[9490]	125	!! = 31 = F(i,j,k,t) = F(local velocity) ( 1/2 \|u\|e laplacian operator
	126	!! or 1/12 \|u\|e^3 bilaplacian operator )
[5836]	127	!! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init
	128	!!
[9490]	129	!! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true
	130	!! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true
[3]	131	!!----------------------------------------------------------------------
[9490]	132	INTEGER :: jk ! dummy loop indices
	133	INTEGER :: ioptio, ierr, inum, ios, inn ! local integer
	134	REAL(wp) :: zah_max, zUfac ! - -
	135	CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s)
[9019]	136	!!
[9526]	137	NAMELIST/namtra_ldf/ ln_traldf_OFF, ln_traldf_lap , ln_traldf_blp , & ! type of operator
	138	& ln_traldf_lev, ln_traldf_hor , ln_traldf_triad, & ! acting direction of the operator
	139	& ln_traldf_iso, ln_traldf_msc , rn_slpmax , & ! option for iso-neutral operator
	140	& ln_triad_iso , ln_botmix_triad, rn_sw_triad , & ! option for triad operator
	141	& nn_aht_ijk_t , rn_Ud , rn_Ld ! lateral eddy coefficient
[3]	142	!!----------------------------------------------------------------------
[5836]	143	!
[9169]	144	IF(lwp) THEN ! control print
	145	WRITE(numout,*)
[9490]	146	WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion'
[9169]	147	WRITE(numout,*) '~~~~~~~~~~~~ '
	148	ENDIF
[9490]	149
[9169]	150	!
[5836]	151	! Choice of lateral tracer physics
	152	! =================================
	153	!
[12276]	154	REWIND( numnam_ref )
[4147]	155	READ ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901)
[11536]	156	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist' )
[12276]	157
	158	REWIND( numnam_cfg )
[4147]	159	READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 )
[11536]	160	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' )
[9169]	161	IF(lwm) WRITE( numond, namtra_ldf )
[5836]	162	!
[1601]	163	IF(lwp) THEN ! control print
[9169]	164	WRITE(numout,*) ' Namelist : namtra_ldf --- lateral mixing parameters (type, direction, coefficients)'
[5836]	165	WRITE(numout,*) ' type :'
[9526]	166	WRITE(numout,*) ' no explicit diffusion ln_traldf_OFF = ', ln_traldf_OFF
[5836]	167	WRITE(numout,*) ' laplacian operator ln_traldf_lap = ', ln_traldf_lap
	168	WRITE(numout,*) ' bilaplacian operator ln_traldf_blp = ', ln_traldf_blp
	169	WRITE(numout,*) ' direction of action :'
	170	WRITE(numout,*) ' iso-level ln_traldf_lev = ', ln_traldf_lev
	171	WRITE(numout,*) ' horizontal (geopotential) ln_traldf_hor = ', ln_traldf_hor
	172	WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso
	173	WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad
[9490]	174	WRITE(numout,*) ' use the Method of Stab. Correction ln_traldf_msc = ', ln_traldf_msc
[5836]	175	WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax
	176	WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso
	177	WRITE(numout,*) ' switching triad or not rn_sw_triad = ', rn_sw_triad
	178	WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad
	179	WRITE(numout,*) ' coefficients :'
	180	WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t
[9490]	181	WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ud = ', rn_Ud, ' m/s'
	182	WRITE(numout,*) ' lateral diffusive length (if cst) rn_Ld = ', rn_Ld, ' m'
[3]	183	ENDIF
[5836]	184	!
	185	!
[9490]	186	! Operator and its acting direction (set nldf_tra)
	187	! =================================
	188	!
	189	nldf_tra = np_ERROR
	190	ioptio = 0
[9526]	191	IF( ln_traldf_OFF ) THEN ; nldf_tra = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF
	192	IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF
	193	IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF
	194	IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
[9490]	195	!
[9526]	196	IF( .NOT.ln_traldf_OFF ) THEN !== direction ==>> type of operator ==!
[9490]	197	ioptio = 0
[9737]	198	IF( ln_traldf_lev ) ioptio = ioptio + 1
	199	IF( ln_traldf_hor ) ioptio = ioptio + 1
	200	IF( ln_traldf_iso ) ioptio = ioptio + 1
	201	IF( ln_traldf_triad ) ioptio = ioptio + 1
	202	IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso/triad)' )
[9490]	203	!
	204	! ! defined the type of lateral diffusion from ln_traldf_... logicals
	205	ierr = 0
	206	IF ( ln_traldf_lap ) THEN ! laplacian operator
	207	IF ( ln_zco ) THEN ! z-coordinate
	208	IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level = horizontal (no rotation)
	209	IF ( ln_traldf_hor ) nldf_tra = np_lap ! iso-level = horizontal (no rotation)
	210	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation)
	211	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation)
	212	ENDIF
	213	IF ( ln_zps ) THEN ! z-coordinate with partial step
	214	IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed
	215	IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation)
	216	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation)
	217	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad (rotation)
	218	ENDIF
	219	IF ( ln_sco ) THEN ! s-coordinate
	220	IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level (no rotation)
	221	IF ( ln_traldf_hor ) nldf_tra = np_lap_i ! horizontal ( rotation)
	222	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation)
	223	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation)
	224	ENDIF
	225	ENDIF
	226	!
	227	IF( ln_traldf_blp ) THEN ! bilaplacian operator
	228	IF ( ln_zco ) THEN ! z-coordinate
	229	IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level = horizontal (no rotation)
	230	IF ( ln_traldf_hor ) nldf_tra = np_blp ! iso-level = horizontal (no rotation)
	231	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
	232	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
	233	ENDIF
	234	IF ( ln_zps ) THEN ! z-coordinate with partial step
	235	IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed
	236	IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation)
	237	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
	238	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
	239	ENDIF
	240	IF ( ln_sco ) THEN ! s-coordinate
	241	IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level (no rotation)
	242	IF ( ln_traldf_hor ) nldf_tra = np_blp_it ! horizontal ( rotation)
	243	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
	244	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
	245	ENDIF
	246	ENDIF
	247	IF ( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' )
[5836]	248	ENDIF
	249	!
[9490]	250	IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) &
	251	& CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' )
	252	IF( ln_isfcav .AND. ln_traldf_triad ) &
	253	& CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' )
	254	!
	255	IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. &
	256	& nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required
	257	!
[5836]	258	IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC
	259	IF( .NOT.ln_traldf_msc ) CALL ctl_stop( 'tra_ldf_init: iso-neutral bilaplacian requires ln_traldf_msc=.true.' )
	260	ENDIF
	261	!
[9490]	262	IF(lwp) THEN
	263	WRITE(numout,*)
	264	SELECT CASE( nldf_tra )
	265	CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion'
	266	CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator'
	267	CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (standard)'
	268	CASE( np_lap_it ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (triad)'
	269	CASE( np_blp ) ; WRITE(numout,*) ' ==>>> bilaplacian iso-level operator'
	270	CASE( np_blp_i ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (standard)'
	271	CASE( np_blp_it ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (triad)'
	272	END SELECT
	273	WRITE(numout,*)
	274	ENDIF
	275
	276	!
[5836]	277	! Space/time variation of eddy coefficients
	278	! ===========================================
	279	!
[9490]	280	l_ldftra_time = .FALSE. ! no time variation except in case defined below
[5836]	281	!
[9526]	282	IF( ln_traldf_OFF ) THEN !== no explicit diffusive operator ==!
[5836]	283	!
[9490]	284	IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated'
	285	RETURN
[5836]	286	!
[9490]	287	ELSE !== a lateral diffusion operator is used ==!
	288	!
	289	! ! allocate the aht arrays
	290	ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr )
	291	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays')
	292	!
	293	ahtu(:,:,jpk) = 0._wp ! last level always 0
	294	ahtv(:,:,jpk) = 0._wp
	295	!.
	296	! ! value of lap/blp eddy mixing coef.
	297	IF( ln_traldf_lap ) THEN ; zUfac = r1_2 *rn_Ud ; inn = 1 ; cl_Units = ' m2/s' ! laplacian
	298	ELSEIF( ln_traldf_blp ) THEN ; zUfac = r1_12*rn_Ud ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian
	299	ENDIF
	300	aht0 = zUfac * rn_Ld**inn ! mixing coefficient
	301	zah_max = zUfac * (rarad)*inn ! maximum reachable coefficient (value at the Equator for e1=1 degree)
	302	!
	303	!
	304	SELECT CASE( nn_aht_ijk_t ) !* Specification of space-time variations of ahtu, ahtv
	305	!
[5836]	306	CASE( 0 ) !== constant ==!
[9490]	307	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = constant = ', aht0, cl_Units
	308	ahtu(:,:,1:jpkm1) = aht0
	309	ahtv(:,:,1:jpkm1) = aht0
[5836]	310	!
	311	CASE( 10 ) !== fixed profile ==!
[9490]	312	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( depth )'
	313	IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, cl_Units
	314	ahtu(:,:,1) = aht0 ! constant surface value
	315	ahtv(:,:,1) = aht0
	316	CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
[5836]	317	!
[9490]	318	CASE ( -20 ) !== fixed horizontal shape and magnitude read in file ==!
	319	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file'
[5836]	320	CALL iom_open( 'eddy_diffusivity_2D.nc', inum )
	321	CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) )
	322	CALL iom_get ( inum, jpdom_data, 'ahtv_2D', ahtv(:,:,1) )
	323	CALL iom_close( inum )
	324	DO jk = 2, jpkm1
[9490]	325	ahtu(:,:,jk) = ahtu(:,:,1)
	326	ahtv(:,:,jk) = ahtv(:,:,1)
[5836]	327	END DO
	328	!
	329	CASE( 20 ) !== fixed horizontal shape ==!
[9490]	330	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)'
	331	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)'
	332	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)'
	333	CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! value proportional to scale factor^inn
[5836]	334	!
	335	CASE( 21 ) !== time varying 2D field ==!
[9490]	336	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, time )'
	337	IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
	338	IF(lwp) WRITE(numout,) ' min value = 0.2 aht0 (with aht0= 1/2 rn_Ud*rn_Ld)'
	339	IF(lwp) WRITE(numout,) ' max value = aei0 (with aei0=1/2 rn_UeLe increased to aht0 within 20N-20S'
[5836]	340	!
	341	l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
	342	!
[9490]	343	IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)', &
	344	& ' incompatible with bilaplacian operator' )
[5836]	345	!
	346	CASE( -30 ) !== fixed 3D shape read in file ==!
[9490]	347	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file'
[5836]	348	CALL iom_open( 'eddy_diffusivity_3D.nc', inum )
	349	CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu )
	350	CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv )
	351	CALL iom_close( inum )
	352	!
	353	CASE( 30 ) !== fixed 3D shape ==!
[9490]	354	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth )'
	355	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)'
	356	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)'
	357	CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! surface value proportional to scale factor^inn
	358	CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) ! reduction with depth
[5836]	359	!
	360	CASE( 31 ) !== time varying 3D field ==!
[9490]	361	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth , time )'
	362	IF(lwp) WRITE(numout,*) ' proportional to the velocity : 1/2 \|u\|e or 1/12 \|u\|e^3'
[5836]	363	!
	364	l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
	365	!
	366	CASE DEFAULT
	367	CALL ctl_stop('ldf_tra_init: wrong choice for nn_aht_ijk_t, the type of space-time variation of aht')
	368	END SELECT
	369	!
[9490]	370	IF( .NOT.l_ldftra_time ) THEN !* No time variation
	371	IF( ln_traldf_lap ) THEN ! laplacian operator (mask only)
	372	ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1)
	373	ahtv(:,:,1:jpkm1) = ahtv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1)
	374	ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator (square root + mask)
	375	ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1)
	376	ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1)
	377	ENDIF
[5836]	378	ENDIF
	379	!
	380	ENDIF
	381	!
	382	END SUBROUTINE ldf_tra_init
[3]	383
	384
[5836]	385	SUBROUTINE ldf_tra( kt )
	386	!!----------------------------------------------------------------------
	387	!! * ROUTINE ldf_tra *
	388	!!
	389	!! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv)
	390	!!
[9490]	391	!! ** Method : * time varying eddy diffusivity coefficients:
[5836]	392	!!
	393	!! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability)
	394	!! with a reduction to 0 in vicinity of the Equator
	395	!! nn_aht_ijk_t = 21 ahtu, ahtv = F(i,j, t) = F(growth rate of baroclinic instability)
	396	!!
	397	!! = 31 ahtu, ahtv = F(i,j,k,t) = F(local velocity) ( \|u\|e /12 laplacian operator
	398	!! or \|u\|e^3/12 bilaplacian operator )
	399	!!
[9490]	400	!! * time varying EIV coefficients: call to ldf_eiv routine
	401	!!
[5836]	402	!! ** action : ahtu, ahtv update at each time step
	403	!! aeiu, aeiv - - - - (if ln_ldfeiv=T)
	404	!!----------------------------------------------------------------------
	405	INTEGER, INTENT(in) :: kt ! time step
	406	!
	407	INTEGER :: ji, jj, jk ! dummy loop indices
[9490]	408	REAL(wp) :: zaht, zahf, zaht_min, zDaht, z1_f20 ! local scalar
[5836]	409	!!----------------------------------------------------------------------
	410	!
[7646]	411	IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients
[5836]	412	! ! =F(growth rate of baroclinic instability)
[9490]	413	! ! max value aeiv_0 ; decreased to 0 within 20N-20S
	414	CALL ldf_eiv( kt, aei0, aeiu, aeiv )
[5836]	415	ENDIF
	416	!
	417	SELECT CASE( nn_aht_ijk_t ) ! Eddy diffusivity coefficients
	418	!
	419	CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability )
[9490]	420	! ! min value 0.2*aht0
	421	! ! max value aht0 (aei0 if nn_aei_ijk_t=21)
	422	! ! increase to aht0 within 20N-20S
[7646]	423	IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei.
[7753]	424	ahtu(:,:,1) = aeiu(:,:,1)
	425	ahtv(:,:,1) = aeiv(:,:,1)
[7646]	426	ELSE ! compute aht.
[9490]	427	CALL ldf_eiv( kt, aht0, ahtu, ahtv )
[5836]	428	ENDIF
	429	!
[9490]	430	z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees)
	431	zaht_min = 0.2_wp * aht0 ! minimum value for aht
	432	zDaht = aht0 - zaht_min
[5836]	433	DO jj = 1, jpj
	434	DO ji = 1, jpi
[7646]	435	!!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg)
	436	!! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points
[9490]	437	zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht
	438	zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht
	439	ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min
	440	ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N
[5836]	441	END DO
	442	END DO
[9490]	443	DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels
[7753]	444	ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)
	445	ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)
[5836]	446	END DO
	447	!
	448	CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
	449	IF( ln_traldf_lap ) THEN ! laplacian operator \|u\| e /12
	450	DO jk = 1, jpkm1
[9490]	451	ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vb are masked
[7753]	452	ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12
[5836]	453	END DO
	454	ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( \|u\| e^3 /12 ) = sqrt( \|u\| e /12 ) * e
	455	DO jk = 1, jpkm1
[7753]	456	ahtu(:,:,jk) = SQRT( ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ) * e1u(:,:)
	457	ahtv(:,:,jk) = SQRT( ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 ) * e2v(:,:)
[5836]	458	END DO
	459	ENDIF
	460	!
	461	END SELECT
	462	!
[7646]	463	CALL iom_put( "ahtu_2d", ahtu(:,:,1) ) ! surface u-eddy diffusivity coeff.
	464	CALL iom_put( "ahtv_2d", ahtv(:,:,1) ) ! surface v-eddy diffusivity coeff.
	465	CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff.
	466	CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff.
	467	!
	468	IF( ln_ldfeiv ) THEN
	469	CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff.
	470	CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff.
	471	CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff.
	472	CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff.
[5836]	473	ENDIF
	474	!
	475	END SUBROUTINE ldf_tra
[1601]	476
[3]	477
[5836]	478	SUBROUTINE ldf_eiv_init
	479	!!----------------------------------------------------------------------
	480	!! * ROUTINE ldf_eiv_init *
	481	!!
	482	!! ** Purpose : initialization of the eiv coeff. from namelist choices.
	483	!!
[9490]	484	!! ** Method : the eiv diffusivity coef. specification depends on:
	485	!! nn_aei_ijk_t = 0 => = constant
	486	!! !
	487	!! = 10 => = F(z) : constant with a reduction of 1/4 with depth
	488	!! !
	489	!! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file
	490	!! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
	491	!! = 21 = F(i,j,t) = F(growth rate of baroclinic instability)
	492	!! !
	493	!! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file
	494	!! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10)
[5836]	495	!!
[9490]	496	!! ** Action : aeiu , aeiv : initialized one for all or l_ldftra_time set to true
	497	!! l_ldfeiv_time : =T if EIV coefficients vary with time
[5836]	498	!!----------------------------------------------------------------------
[9490]	499	INTEGER :: jk ! dummy loop indices
	500	INTEGER :: ierr, inum, ios, inn ! local integer
	501	REAL(wp) :: zah_max, zUfac ! - scalar
	502	!!
	503	NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv)
	504	& nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient
[5836]	505	!!----------------------------------------------------------------------
	506	!
[9169]	507	IF(lwp) THEN ! control print
	508	WRITE(numout,*)
	509	WRITE(numout,*) 'ldf_eiv_init : eddy induced velocity parametrization'
	510	WRITE(numout,*) '~~~~~~~~~~~~ '
	511	ENDIF
	512	!
[9490]	513	READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901)
[11536]	514	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' )
[5836]	515	!
[9490]	516	READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 )
[11536]	517	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' )
[9490]	518	IF(lwm) WRITE ( numond, namtra_eiv )
[5836]	519
	520	IF(lwp) THEN ! control print
[9490]	521	WRITE(numout,*) ' Namelist namtra_eiv : '
	522	WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv
	523	WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia
	524	WRITE(numout,*) ' coefficients :'
	525	WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aht_ijk_t
	526	WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s'
	527	WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m'
[5836]	528	WRITE(numout,*)
[3]	529	ENDIF
[5836]	530	!
[9490]	531	l_ldfeiv_time = .FALSE. ! no time variation except in case defined below
[5836]	532	!
[9490]	533	!
	534	IF( .NOT.ln_ldfeiv ) THEN !== Parametrization not used ==!
	535	!
	536	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used'
	537	ln_ldfeiv_dia = .FALSE.
	538	!
	539	ELSE !== use the parametrization ==!
	540	!
	541	IF(lwp) WRITE(numout,*) ' ==>>> use eddy induced velocity parametrization'
	542	IF(lwp) WRITE(numout,*)
	543	!
	544	IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' )
	545	!
	546	! != allocate the aei arrays
[5836]	547	ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr )
	548	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays')
	549	!
[9490]	550	! != Specification of space-time variations of eaiu, aeiv
[5836]	551	!
[9490]	552	aeiu(:,:,jpk) = 0._wp ! last level always 0
	553	aeiv(:,:,jpk) = 0._wp
	554	! ! value of EIV coef. (laplacian operator)
	555	zUfac = r1_2 *rn_Ue ! velocity factor
	556	inn = 1 ! L-exponent
	557	aei0 = zUfac * rn_Le**inn ! mixing coefficient
	558	zah_max = zUfac * (rarad)*inn ! maximum reachable coefficient (value at the Equator)
	559
	560	SELECT CASE( nn_aei_ijk_t ) !* Specification of space-time variations
	561	!
	562	CASE( 0 ) !-- constant --!
	563	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', aei0, ' m2/s'
	564	aeiu(:,:,1:jpkm1) = aei0
	565	aeiv(:,:,1:jpkm1) = aei0
[5836]	566	!
[9490]	567	CASE( 10 ) !-- fixed profile --!
[9169]	568	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( depth )'
[9490]	569	IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, ' m2/s'
	570	aeiu(:,:,1) = aei0 ! constant surface value
	571	aeiv(:,:,1) = aei0
	572	CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
[5836]	573	!
[9490]	574	CASE ( -20 ) !-- fixed horizontal shape read in file --!
	575	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file'
[5836]	576	CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum )
	577	CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) )
	578	CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) )
	579	CALL iom_close( inum )
[9490]	580	DO jk = 2, jpkm1
[7753]	581	aeiu(:,:,jk) = aeiu(:,:,1)
	582	aeiv(:,:,jk) = aeiv(:,:,1)
[5836]	583	END DO
	584	!
[9490]	585	CASE( 20 ) !-- fixed horizontal shape --!
	586	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( e1, e2 )'
	587	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ue, ' m/s and Le = Max(e1,e2)'
	588	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s for e1=1°)'
	589	CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! value proportional to scale factor^inn
[5836]	590	!
[9490]	591	CASE( 21 ) !-- time varying 2D field --!
	592	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )'
	593	IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
	594	IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s'
[5836]	595	!
	596	l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
	597	!
[9490]	598	CASE( -30 ) !-- fixed 3D shape read in file --!
	599	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
[5836]	600	CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum )
	601	CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu )
	602	CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv )
	603	CALL iom_close( inum )
	604	!
[9490]	605	CASE( 30 ) !-- fixed 3D shape --!
	606	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, depth )'
	607	CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! surface value proportional to scale factor^inn
	608	CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) ! reduction with depth
[5836]	609	!
	610	CASE DEFAULT
	611	CALL ctl_stop('ldf_tra_init: wrong choice for nn_aei_ijk_t, the type of space-time variation of aei')
	612	END SELECT
	613	!
[9490]	614	IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation
	615	DO jk = 1, jpkm1
	616	aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk)
	617	ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk)
	618	END DO
	619	ENDIF
	620	!
[3]	621	ENDIF
[5836]	622	!
	623	END SUBROUTINE ldf_eiv_init
[3]	624
	625
[5836]	626	SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
	627	!!----------------------------------------------------------------------
	628	!! * ROUTINE ldf_eiv *
	629	!!
	630	!! ** Purpose : Compute the eddy induced velocity coefficient from the
	631	!! growth rate of baroclinic instability.
	632	!!
	633	!! ** Method : coefficient function of the growth rate of baroclinic instability
	634	!!
	635	!! Reference : Treguier et al. JPO 1997 ; Held and Larichev JAS 1996
	636	!!----------------------------------------------------------------------
	637	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	638	REAL(wp) , INTENT(inout) :: paei0 ! max value [m2/s]
	639	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s]
	640	!
	641	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	642	REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars
[9490]	643	REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace
[5836]	644	!!----------------------------------------------------------------------
	645	!
[9490]	646	zn (:,:) = 0._wp ! Local initialization
	647	zhw(:,:) = 5._wp
	648	zah(:,:) = 0._wp
	649	zRo(:,:) = 0._wp
[5836]	650	! ! Compute lateral diffusive coefficient at T-point
	651	IF( ln_traldf_triad ) THEN
	652	DO jk = 1, jpk
	653	DO jj = 2, jpjm1
	654	DO ji = 2, jpim1
	655	! Take the max of N^2 and zero then take the vertical sum
	656	! of the square root of the resulting N^2 ( required to compute
	657	! internal Rossby radius Ro = .5 * sum_jpk(N) / f
	658	zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
[6140]	659	zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
[5836]	660	! Compute elements required for the inverse time scale of baroclinic
	661	! eddies using the isopycnal slopes calculated in ldfslp.F :
	662	! T^-1 = sqrt(m_jpk(N^2(r1^2+r2^2)e3w))
[12276]	663	ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
[5836]	664	zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
	665	zhw(ji,jj) = zhw(ji,jj) + ze3w
	666	END DO
	667	END DO
	668	END DO
	669	ELSE
	670	DO jk = 1, jpk
	671	DO jj = 2, jpjm1
	672	DO ji = 2, jpim1
	673	! Take the max of N^2 and zero then take the vertical sum
	674	! of the square root of the resulting N^2 ( required to compute
	675	! internal Rossby radius Ro = .5 * sum_jpk(N) / f
	676	zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
[6140]	677	zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
[5836]	678	! Compute elements required for the inverse time scale of baroclinic
	679	! eddies using the isopycnal slopes calculated in ldfslp.F :
	680	! T^-1 = sqrt(m_jpk(N^2(r1^2+r2^2)e3w))
[12276]	681	ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
[5836]	682	zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	683	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
	684	zhw(ji,jj) = zhw(ji,jj) + ze3w
	685	END DO
	686	END DO
	687	END DO
[9490]	688	ENDIF
[5836]	689
	690	DO jj = 2, jpjm1
	691	DO ji = fs_2, fs_jpim1 ! vector opt.
	692	zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
[9490]	693	! Rossby radius at w-point taken betwenn 2 km and 40km
	694	zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) )
[5836]	695	! Compute aeiw by multiplying Ro^2 and T^-1
[9490]	696	zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
[5836]	697	END DO
	698	END DO
	699
	700	! !== Bound on eiv coeff. ==!
	701	z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) )
	702	DO jj = 2, jpjm1
	703	DO ji = fs_2, fs_jpim1 ! vector opt.
[9490]	704	zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
[5836]	705	zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0
	706	END DO
	707	END DO
[10425]	708	CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. ) ! lateral boundary condition
[5836]	709	!
	710	DO jj = 2, jpjm1 !== aei at u- and v-points ==!
	711	DO ji = fs_2, fs_jpim1 ! vector opt.
	712	paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1)
	713	paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1)
	714	END DO
	715	END DO
[10425]	716	CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition
[5836]	717
	718	DO jk = 2, jpkm1 !== deeper values equal the surface one ==!
[7753]	719	paeiu(:,:,jk) = paeiu(:,:,1) * umask(:,:,jk)
	720	paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk)
[5836]	721	END DO
	722	!
	723	END SUBROUTINE ldf_eiv
	724
	725
	726	SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype )
	727	!!----------------------------------------------------------------------
	728	!! * ROUTINE ldf_eiv_trp *
	729	!!
	730	!! ** Purpose : add to the input ocean transport the contribution of
	731	!! the eddy induced velocity parametrization.
	732	!!
	733	!! ** Method : The eddy induced transport is computed from a flux stream-
	734	!! function which depends on the slope of iso-neutral surfaces
	735	!! (see ldf_slp). For example, in the i-k plan :
	736	!! psi_uw = mk(aeiu) e2u mi(wslpi) [in m3/s]
	737	!! Utr_eiv = - dk[psi_uw]
	738	!! Vtr_eiv = + di[psi_uw]
	739	!! ln_ldfeiv_dia = T : output the associated streamfunction,
	740	!! velocity and heat transport (call ldf_eiv_dia)
	741	!!
	742	!! ** Action : pun, pvn increased by the eiv transport
	743	!!----------------------------------------------------------------------
	744	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	745	INTEGER , INTENT(in ) :: kit000 ! first time step index
	746	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	747	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun ! in : 3 ocean transport components [m3/s]
	748	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pvn ! out: 3 ocean transport components [m3/s]
	749	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pwn ! increased by the eiv [m3/s]
	750	!!
	751	INTEGER :: ji, jj, jk ! dummy loop indices
	752	REAL(wp) :: zuwk, zuwk1, zuwi, zuwi1 ! local scalars
	753	REAL(wp) :: zvwk, zvwk1, zvwj, zvwj1 ! - -
[9019]	754	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
[5836]	755	!!----------------------------------------------------------------------
	756	!
	757	IF( kt == kit000 ) THEN
	758	IF(lwp) WRITE(numout,*)
	759	IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :'
	760	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ add to velocity fields the eiv component'
[3]	761	ENDIF
[3634]	762
[5836]	763
[7753]	764	zpsi_uw(:,:, 1 ) = 0._wp ; zpsi_vw(:,:, 1 ) = 0._wp
	765	zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp
[5836]	766	!
	767	DO jk = 2, jpkm1
	768	DO jj = 1, jpjm1
	769	DO ji = 1, fs_jpim1 ! vector opt.
[9490]	770	zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) &
	771	& * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk)
	772	zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) &
	773	& * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk)
[5836]	774	END DO
	775	END DO
	776	END DO
	777	!
	778	DO jk = 1, jpkm1
	779	DO jj = 1, jpjm1
	780	DO ji = 1, fs_jpim1 ! vector opt.
	781	pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
	782	pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
	783	END DO
	784	END DO
	785	END DO
	786	DO jk = 1, jpkm1
	787	DO jj = 2, jpjm1
	788	DO ji = fs_2, fs_jpim1 ! vector opt.
	789	pwn(ji,jj,jk) = pwn(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) &
	790	& + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) )
	791	END DO
	792	END DO
	793	END DO
	794	!
	795	! ! diagnose the eddy induced velocity and associated heat transport
	796	IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
	797	!
	798	END SUBROUTINE ldf_eiv_trp
[3634]	799
[5836]	800
	801	SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )
	802	!!----------------------------------------------------------------------
	803	!! * ROUTINE ldf_eiv_dia *
	804	!!
	805	!! ** Purpose : diagnose the eddy induced velocity and its associated
	806	!! vertically integrated heat transport.
	807	!!
	808	!! ** Method :
	809	!!
	810	!!----------------------------------------------------------------------
	811	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s]
[1601]	812	!
[5836]	813	INTEGER :: ji, jj, jk ! dummy loop indices
	814	REAL(wp) :: zztmp ! local scalar
[9019]	815	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	816	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[5836]	817	!!----------------------------------------------------------------------
	818	!
[9019]	819	!!gm I don't like this routine.... Crazy way of doing things, not optimal at all...
	820	!!gm to be redesigned....
[5836]	821	! !== eiv stream function: output ==!
[10425]	822	CALL lbc_lnk_multi( 'ldftra', psi_uw, 'U', -1. , psi_vw, 'V', -1. )
[5836]	823	!
	824	!!gm CALL iom_put( "psi_eiv_uw", psi_uw ) ! output
	825	!!gm CALL iom_put( "psi_eiv_vw", psi_vw )
	826	!
	827	! !== eiv velocities: calculate and output ==!
	828	!
[7753]	829	zw3d(:,:,jpk) = 0._wp ! bottom value always 0
[5836]	830	!
	831	DO jk = 1, jpkm1 ! e2u e3u u_eiv = -dk[psi_uw]
[7753]	832	zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
[5836]	833	END DO
	834	CALL iom_put( "uoce_eiv", zw3d )
	835	!
	836	DO jk = 1, jpkm1 ! e1v e3v v_eiv = -dk[psi_vw]
[7753]	837	zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
[5836]	838	END DO
	839	CALL iom_put( "voce_eiv", zw3d )
	840	!
	841	DO jk = 1, jpkm1 ! e1 e2 w_eiv = dk[psix] + dk[psix]
	842	DO jj = 2, jpjm1
	843	DO ji = fs_2, fs_jpim1 ! vector opt.
	844	zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) &
	845	& + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj)
	846	END DO
	847	END DO
	848	END DO
[10425]	849	CALL lbc_lnk( 'ldftra', zw3d, 'T', 1. ) ! lateral boundary condition
[5836]	850	CALL iom_put( "woce_eiv", zw3d )
	851	!
[12276]	852	IF( iom_use('weiv_masstr') ) THEN ! vertical mass transport & its square value
	853	zw2d(:,:) = rau0 * e1e2t(:,:)
	854	DO jk = 1, jpk
	855	zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
	856	END DO
	857	CALL iom_put( "weiv_masstr" , zw3d )
	858	ENDIF
[5836]	859	!
[12276]	860	IF( iom_use('ueiv_masstr') ) THEN
	861	zw3d(:,:,:) = 0.e0
	862	DO jk = 1, jpkm1
	863	zw3d(:,:,jk) = rau0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) )
	864	END DO
	865	CALL iom_put( "ueiv_masstr", zw3d ) ! mass transport in i-direction
	866	ENDIF
	867	!
[7646]	868	zztmp = 0.5_wp * rau0 * rcp
	869	IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
[7753]	870	zw2d(:,:) = 0._wp
	871	zw3d(:,:,:) = 0._wp
	872	DO jk = 1, jpkm1
	873	DO jj = 2, jpjm1
	874	DO ji = fs_2, fs_jpim1 ! vector opt.
	875	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) &
	876	& * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) )
	877	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
	878	END DO
	879	END DO
	880	END DO
[10425]	881	CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. )
	882	CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. )
[7753]	883	CALL iom_put( "ueiv_heattr" , zztmp * zw2d ) ! heat transport in i-direction
	884	CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction
[7646]	885	ENDIF
[12276]	886	!
	887	IF( iom_use('veiv_masstr') ) THEN
	888	zw3d(:,:,:) = 0.e0
	889	DO jk = 1, jpkm1
	890	zw3d(:,:,jk) = rau0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) )
	891	END DO
	892	CALL iom_put( "veiv_masstr", zw3d ) ! mass transport in i-direction
	893	ENDIF
	894	!
[7753]	895	zw2d(:,:) = 0._wp
	896	zw3d(:,:,:) = 0._wp
[7646]	897	DO jk = 1, jpkm1
	898	DO jj = 2, jpjm1
	899	DO ji = fs_2, fs_jpim1 ! vector opt.
	900	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) &
	901	& * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) )
	902	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
[5836]	903	END DO
	904	END DO
[7646]	905	END DO
[10425]	906	CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. )
[7646]	907	CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction
	908	CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction
	909	!
[12276]	910	IF( iom_use( 'sophteiv' ) ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
[7646]	911	!
	912	zztmp = 0.5_wp * 0.5
	913	IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
[7753]	914	zw2d(:,:) = 0._wp
	915	zw3d(:,:,:) = 0._wp
	916	DO jk = 1, jpkm1
	917	DO jj = 2, jpjm1
	918	DO ji = fs_2, fs_jpim1 ! vector opt.
	919	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) &
	920	& * ( tsn (ji,jj,jk,jp_sal) + tsn (ji+1,jj,jk,jp_sal) )
	921	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
	922	END DO
	923	END DO
	924	END DO
[10425]	925	CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. )
	926	CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. )
[7753]	927	CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction
	928	CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction
[7646]	929	ENDIF
[7753]	930	zw2d(:,:) = 0._wp
	931	zw3d(:,:,:) = 0._wp
[7646]	932	DO jk = 1, jpkm1
	933	DO jj = 2, jpjm1
	934	DO ji = fs_2, fs_jpim1 ! vector opt.
	935	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) &
	936	& * ( tsn (ji,jj,jk,jp_sal) + tsn (ji,jj+1,jk,jp_sal) )
	937	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
[5836]	938	END DO
	939	END DO
[7646]	940	END DO
[10425]	941	CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. )
[7646]	942	CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction
	943	CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction
[5836]	944	!
[12276]	945	IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
[7646]	946	!
	947	!
[5836]	948	END SUBROUTINE ldf_eiv_dia
[3]	949
	950	!!======================================================================
	951	END MODULE ldftra

Note: See TracBrowser for help on using the repository browser.

Download in other formats: