Context Navigation

ldftra.F90 @ 9526

Last change on this file since 9526 was 9526, checked in by gm, 6 years ago
dev_merge_2017: rename ln_...NONE as ln_...OFF (CONFIG, OPA_SRC, TOP_SRC) ; agrif zoom now only in AGRIF_NORDIC configuration
Property svn:keywords set to `Id`
File size: 52.5 KB

Line
1	MODULE ldftra
2	!!======================================================================
3	!! * MODULE ldftra *
4	!! Ocean physics: lateral diffusivity coefficients
5	!!=====================================================================
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
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
14	!! ldf_tra_init : initialization, namelist read, and parameters control
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
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers
22	USE dom_oce ! ocean space and time domain
23	USE phycst ! physical constants
24	USE ldfslp ! lateral diffusion: slope of iso-neutral surfaces
25	USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases
26	USE diaptr
27	!
28	USE in_out_manager ! I/O manager
29	USE iom ! I/O module for ehanced bottom friction file
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
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 =!
45	LOGICAL , PUBLIC :: ln_traldf_OFF !: no operator: No explicit diffusion
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)
52	! != iso-neutral options =!
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.
61	! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case)
62	! ! bht_0 = 1/12 Ud*Ld^3 (blp case)
63	REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s]
64	REAL(wp), PUBLIC :: rn_Ld !: lateral diffusive length [m]
65
66	! !!* Namelist namtra_eiv : eddy induced velocity param. *
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.
72	REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s]
73	REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m]
74
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)
84	LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef.
85	LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef.
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
90	REAL(wp) :: aht0, aei0 ! constant eddy coefficients (deduced from namelist values) [m2/s]
91	REAL(wp) :: r1_2 = 0.5_wp ! =1/2
92	REAL(wp) :: r1_4 = 0.25_wp ! =1/4
93	REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12
94
95	!! * Substitutions
96	# include "vectopt_loop_substitute.h90"
97	!!----------------------------------------------------------------------
98	!! NEMO/OPA 3.7 , NEMO Consortium (2015)
99	!! $Id$
100	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
101	!!----------------------------------------------------------------------
102	CONTAINS
103
104	SUBROUTINE ldf_tra_init
105	!!----------------------------------------------------------------------
106	!! * ROUTINE ldf_tra_init *
107	!!
108	!! ** Purpose : initializations of the tracer lateral mixing coeff.
109	!!
110	!! ** Method : * the eddy diffusivity coef. specification depends on:
111	!!
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)
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 )
127	!! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init
128	!!
129	!! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true
130	!! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true
131	!!----------------------------------------------------------------------
132	INTEGER :: jk ! dummy loop indices
133	INTEGER :: ioptio, ierr, inum, ios, inn ! local integer
134	REAL(wp) :: zah_max, zUfac ! - -
135	CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s)
136	!!
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
142	!!----------------------------------------------------------------------
143	!
144	IF(lwp) THEN ! control print
145	WRITE(numout,*)
146	WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion'
147	WRITE(numout,*) '~~~~~~~~~~~~ '
148	ENDIF
149
150	!
151	! Choice of lateral tracer physics
152	! =================================
153	!
154	REWIND( numnam_ref ) ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers
155	READ ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901)
156	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist', lwp )
157	REWIND( numnam_cfg ) ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers
158	READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 )
159	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist', lwp )
160	IF(lwm) WRITE( numond, namtra_ldf )
161	!
162	IF(lwp) THEN ! control print
163	WRITE(numout,*) ' Namelist : namtra_ldf --- lateral mixing parameters (type, direction, coefficients)'
164	WRITE(numout,*) ' type :'
165	WRITE(numout,*) ' no explicit diffusion ln_traldf_OFF = ', ln_traldf_OFF
166	WRITE(numout,*) ' laplacian operator ln_traldf_lap = ', ln_traldf_lap
167	WRITE(numout,*) ' bilaplacian operator ln_traldf_blp = ', ln_traldf_blp
168	WRITE(numout,*) ' direction of action :'
169	WRITE(numout,*) ' iso-level ln_traldf_lev = ', ln_traldf_lev
170	WRITE(numout,*) ' horizontal (geopotential) ln_traldf_hor = ', ln_traldf_hor
171	WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso
172	WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad
173	WRITE(numout,*) ' use the Method of Stab. Correction ln_traldf_msc = ', ln_traldf_msc
174	WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax
175	WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso
176	WRITE(numout,*) ' switching triad or not rn_sw_triad = ', rn_sw_triad
177	WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad
178	WRITE(numout,*) ' coefficients :'
179	WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t
180	WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ud = ', rn_Ud, ' m/s'
181	WRITE(numout,*) ' lateral diffusive length (if cst) rn_Ld = ', rn_Ld, ' m'
182	ENDIF
183	!
184	!
185	! Operator and its acting direction (set nldf_tra)
186	! =================================
187	!
188	nldf_tra = np_ERROR
189	ioptio = 0
190	IF( ln_traldf_OFF ) THEN ; nldf_tra = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF
191	IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF
192	IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF
193	IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
194	!
195	IF( .NOT.ln_traldf_OFF ) THEN !== direction ==>> type of operator ==!
196	ioptio = 0
197	IF( ln_traldf_lev ) ioptio = ioptio + 1
198	IF( ln_traldf_hor ) ioptio = ioptio + 1
199	IF( ln_traldf_iso ) ioptio = ioptio + 1
200	IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' )
201	!
202	! ! defined the type of lateral diffusion from ln_traldf_... logicals
203	ierr = 0
204	IF ( ln_traldf_lap ) THEN ! laplacian operator
205	IF ( ln_zco ) THEN ! z-coordinate
206	IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level = horizontal (no rotation)
207	IF ( ln_traldf_hor ) nldf_tra = np_lap ! iso-level = horizontal (no rotation)
208	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation)
209	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation)
210	ENDIF
211	IF ( ln_zps ) THEN ! z-coordinate with partial step
212	IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed
213	IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation)
214	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation)
215	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad (rotation)
216	ENDIF
217	IF ( ln_sco ) THEN ! s-coordinate
218	IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level (no rotation)
219	IF ( ln_traldf_hor ) nldf_tra = np_lap_i ! horizontal ( rotation)
220	IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation)
221	IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation)
222	ENDIF
223	ENDIF
224	!
225	IF( ln_traldf_blp ) THEN ! bilaplacian operator
226	IF ( ln_zco ) THEN ! z-coordinate
227	IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level = horizontal (no rotation)
228	IF ( ln_traldf_hor ) nldf_tra = np_blp ! iso-level = horizontal (no rotation)
229	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
230	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
231	ENDIF
232	IF ( ln_zps ) THEN ! z-coordinate with partial step
233	IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed
234	IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation)
235	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
236	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
237	ENDIF
238	IF ( ln_sco ) THEN ! s-coordinate
239	IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level (no rotation)
240	IF ( ln_traldf_hor ) nldf_tra = np_blp_it ! horizontal ( rotation)
241	IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation)
242	IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation)
243	ENDIF
244	ENDIF
245	IF ( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' )
246	ENDIF
247	!
248	IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) &
249	& CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' )
250	IF( ln_isfcav .AND. ln_traldf_triad ) &
251	& CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' )
252	!
253	IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. &
254	& nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required
255	!
256	IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC
257	IF( .NOT.ln_traldf_msc ) CALL ctl_stop( 'tra_ldf_init: iso-neutral bilaplacian requires ln_traldf_msc=.true.' )
258	ENDIF
259	!
260	IF(lwp) THEN
261	WRITE(numout,*)
262	SELECT CASE( nldf_tra )
263	CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion'
264	CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator'
265	CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (standard)'
266	CASE( np_lap_it ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (triad)'
267	CASE( np_blp ) ; WRITE(numout,*) ' ==>>> bilaplacian iso-level operator'
268	CASE( np_blp_i ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (standard)'
269	CASE( np_blp_it ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (triad)'
270	END SELECT
271	WRITE(numout,*)
272	ENDIF
273
274	!
275	! Space/time variation of eddy coefficients
276	! ===========================================
277	!
278	l_ldftra_time = .FALSE. ! no time variation except in case defined below
279	!
280	IF( ln_traldf_OFF ) THEN !== no explicit diffusive operator ==!
281	!
282	IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated'
283	RETURN
284	!
285	ELSE !== a lateral diffusion operator is used ==!
286	!
287	! ! allocate the aht arrays
288	ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr )
289	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays')
290	!
291	ahtu(:,:,jpk) = 0._wp ! last level always 0
292	ahtv(:,:,jpk) = 0._wp
293	!.
294	! ! value of lap/blp eddy mixing coef.
295	IF( ln_traldf_lap ) THEN ; zUfac = r1_2 *rn_Ud ; inn = 1 ; cl_Units = ' m2/s' ! laplacian
296	ELSEIF( ln_traldf_blp ) THEN ; zUfac = r1_12*rn_Ud ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian
297	ENDIF
298	aht0 = zUfac * rn_Ld**inn ! mixing coefficient
299	zah_max = zUfac * (rarad)*inn ! maximum reachable coefficient (value at the Equator for e1=1 degree)
300	!
301	!
302	SELECT CASE( nn_aht_ijk_t ) !* Specification of space-time variations of ahtu, ahtv
303	!
304	CASE( 0 ) !== constant ==!
305	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = constant = ', aht0, cl_Units
306	ahtu(:,:,1:jpkm1) = aht0
307	ahtv(:,:,1:jpkm1) = aht0
308	!
309	CASE( 10 ) !== fixed profile ==!
310	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( depth )'
311	IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, cl_Units
312	ahtu(:,:,1) = aht0 ! constant surface value
313	ahtv(:,:,1) = aht0
314	CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
315	!
316	CASE ( -20 ) !== fixed horizontal shape and magnitude read in file ==!
317	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file'
318	CALL iom_open( 'eddy_diffusivity_2D.nc', inum )
319	CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) )
320	CALL iom_get ( inum, jpdom_data, 'ahtv_2D', ahtv(:,:,1) )
321	CALL iom_close( inum )
322	DO jk = 2, jpkm1
323	ahtu(:,:,jk) = ahtu(:,:,1)
324	ahtv(:,:,jk) = ahtv(:,:,1)
325	END DO
326	!
327	CASE( 20 ) !== fixed horizontal shape ==!
328	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)'
329	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)'
330	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)'
331	CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! value proportional to scale factor^inn
332	!
333	CASE( 21 ) !== time varying 2D field ==!
334	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, time )'
335	IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
336	IF(lwp) WRITE(numout,) ' min value = 0.2 aht0 (with aht0= 1/2 rn_Ud*rn_Ld)'
337	IF(lwp) WRITE(numout,) ' max value = aei0 (with aei0=1/2 rn_UeLe increased to aht0 within 20N-20S'
338	!
339	l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
340	!
341	IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)', &
342	& ' incompatible with bilaplacian operator' )
343	!
344	CASE( -30 ) !== fixed 3D shape read in file ==!
345	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file'
346	CALL iom_open( 'eddy_diffusivity_3D.nc', inum )
347	CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu )
348	CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv )
349	CALL iom_close( inum )
350	!
351	CASE( 30 ) !== fixed 3D shape ==!
352	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth )'
353	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)'
354	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)'
355	CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! surface value proportional to scale factor^inn
356	CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) ! reduction with depth
357	!
358	CASE( 31 ) !== time varying 3D field ==!
359	IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth , time )'
360	IF(lwp) WRITE(numout,*) ' proportional to the velocity : 1/2 \|u\|e or 1/12 \|u\|e^3'
361	!
362	l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
363	!
364	CASE DEFAULT
365	CALL ctl_stop('ldf_tra_init: wrong choice for nn_aht_ijk_t, the type of space-time variation of aht')
366	END SELECT
367	!
368	IF( .NOT.l_ldftra_time ) THEN !* No time variation
369	IF( ln_traldf_lap ) THEN ! laplacian operator (mask only)
370	ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1)
371	ahtv(:,:,1:jpkm1) = ahtv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1)
372	ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator (square root + mask)
373	ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1)
374	ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1)
375	ENDIF
376	ENDIF
377	!
378	ENDIF
379	!
380	END SUBROUTINE ldf_tra_init
381
382
383	SUBROUTINE ldf_tra( kt )
384	!!----------------------------------------------------------------------
385	!! * ROUTINE ldf_tra *
386	!!
387	!! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv)
388	!!
389	!! ** Method : * time varying eddy diffusivity coefficients:
390	!!
391	!! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability)
392	!! with a reduction to 0 in vicinity of the Equator
393	!! nn_aht_ijk_t = 21 ahtu, ahtv = F(i,j, t) = F(growth rate of baroclinic instability)
394	!!
395	!! = 31 ahtu, ahtv = F(i,j,k,t) = F(local velocity) ( \|u\|e /12 laplacian operator
396	!! or \|u\|e^3/12 bilaplacian operator )
397	!!
398	!! * time varying EIV coefficients: call to ldf_eiv routine
399	!!
400	!! ** action : ahtu, ahtv update at each time step
401	!! aeiu, aeiv - - - - (if ln_ldfeiv=T)
402	!!----------------------------------------------------------------------
403	INTEGER, INTENT(in) :: kt ! time step
404	!
405	INTEGER :: ji, jj, jk ! dummy loop indices
406	REAL(wp) :: zaht, zahf, zaht_min, zDaht, z1_f20 ! local scalar
407	!!----------------------------------------------------------------------
408	!
409	IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients
410	! ! =F(growth rate of baroclinic instability)
411	! ! max value aeiv_0 ; decreased to 0 within 20N-20S
412	CALL ldf_eiv( kt, aei0, aeiu, aeiv )
413	ENDIF
414	!
415	SELECT CASE( nn_aht_ijk_t ) ! Eddy diffusivity coefficients
416	!
417	CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability )
418	! ! min value 0.2*aht0
419	! ! max value aht0 (aei0 if nn_aei_ijk_t=21)
420	! ! increase to aht0 within 20N-20S
421	IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei.
422	ahtu(:,:,1) = aeiu(:,:,1)
423	ahtv(:,:,1) = aeiv(:,:,1)
424	ELSE ! compute aht.
425	CALL ldf_eiv( kt, aht0, ahtu, ahtv )
426	ENDIF
427	!
428	z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees)
429	zaht_min = 0.2_wp * aht0 ! minimum value for aht
430	zDaht = aht0 - zaht_min
431	DO jj = 1, jpj
432	DO ji = 1, jpi
433	!!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg)
434	!! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points
435	zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht
436	zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht
437	ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min
438	ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N
439	END DO
440	END DO
441	DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels
442	ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)
443	ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)
444	END DO
445	!
446	CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
447	IF( ln_traldf_lap ) THEN ! laplacian operator \|u\| e /12
448	DO jk = 1, jpkm1
449	ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vb are masked
450	ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12
451	END DO
452	ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( \|u\| e^3 /12 ) = sqrt( \|u\| e /12 ) * e
453	DO jk = 1, jpkm1
454	ahtu(:,:,jk) = SQRT( ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ) * e1u(:,:)
455	ahtv(:,:,jk) = SQRT( ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 ) * e2v(:,:)
456	END DO
457	ENDIF
458	!
459	END SELECT
460	!
461	CALL iom_put( "ahtu_2d", ahtu(:,:,1) ) ! surface u-eddy diffusivity coeff.
462	CALL iom_put( "ahtv_2d", ahtv(:,:,1) ) ! surface v-eddy diffusivity coeff.
463	CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff.
464	CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff.
465	!
466	IF( ln_ldfeiv ) THEN
467	CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff.
468	CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff.
469	CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff.
470	CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff.
471	ENDIF
472	!
473	END SUBROUTINE ldf_tra
474
475
476	SUBROUTINE ldf_eiv_init
477	!!----------------------------------------------------------------------
478	!! * ROUTINE ldf_eiv_init *
479	!!
480	!! ** Purpose : initialization of the eiv coeff. from namelist choices.
481	!!
482	!! ** Method : the eiv diffusivity coef. specification depends on:
483	!! nn_aei_ijk_t = 0 => = constant
484	!! !
485	!! = 10 => = F(z) : constant with a reduction of 1/4 with depth
486	!! !
487	!! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file
488	!! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
489	!! = 21 = F(i,j,t) = F(growth rate of baroclinic instability)
490	!! !
491	!! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file
492	!! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10)
493	!!
494	!! ** Action : aeiu , aeiv : initialized one for all or l_ldftra_time set to true
495	!! l_ldfeiv_time : =T if EIV coefficients vary with time
496	!!----------------------------------------------------------------------
497	INTEGER :: jk ! dummy loop indices
498	INTEGER :: ierr, inum, ios, inn ! local integer
499	REAL(wp) :: zah_max, zUfac ! - scalar
500	!!
501	NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv)
502	& nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient
503	!!----------------------------------------------------------------------
504	!
505	IF(lwp) THEN ! control print
506	WRITE(numout,*)
507	WRITE(numout,*) 'ldf_eiv_init : eddy induced velocity parametrization'
508	WRITE(numout,*) '~~~~~~~~~~~~ '
509	ENDIF
510	!
511	REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param.
512	READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901)
513	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp )
514	!
515	REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param.
516	READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 )
517	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist', lwp )
518	IF(lwm) WRITE ( numond, namtra_eiv )
519
520	IF(lwp) THEN ! control print
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'
528	WRITE(numout,*)
529	ENDIF
530	!
531	l_ldfeiv_time = .FALSE. ! no time variation except in case defined below
532	!
533	!
534	IF( .NOT.ln_ldfeiv ) THEN !== Parametrization not used ==!
535	!
536	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used'
537	ln_ldfeiv_dia = .FALSE.
538	!
539	ELSE !== use the parametrization ==!
540	!
541	IF(lwp) WRITE(numout,*) ' ==>>> use eddy induced velocity parametrization'
542	IF(lwp) WRITE(numout,*)
543	!
544	IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' )
545	!
546	! != allocate the aei arrays
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	!
550	! != Specification of space-time variations of eaiu, aeiv
551	!
552	aeiu(:,:,jpk) = 0._wp ! last level always 0
553	aeiv(:,:,jpk) = 0._wp
554	! ! value of EIV coef. (laplacian operator)
555	zUfac = r1_2 *rn_Ue ! velocity factor
556	inn = 1 ! L-exponent
557	aei0 = zUfac * rn_Le**inn ! mixing coefficient
558	zah_max = zUfac * (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
566	!
567	CASE( 10 ) !-- fixed profile --!
568	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( depth )'
569	IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, ' m2/s'
570	aeiu(:,:,1) = aei0 ! constant surface value
571	aeiv(:,:,1) = aei0
572	CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
573	!
574	CASE ( -20 ) !-- fixed horizontal shape read in file --!
575	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file'
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 )
580	DO jk = 2, jpkm1
581	aeiu(:,:,jk) = aeiu(:,:,1)
582	aeiv(:,:,jk) = aeiv(:,:,1)
583	END DO
584	!
585	CASE( 20 ) !-- fixed horizontal shape --!
586	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( e1, e2 )'
587	IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ue, ' m/s and Le = Max(e1,e2)'
588	IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s for e1=1°)'
589	CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! value proportional to scale factor^inn
590	!
591	CASE( 21 ) !-- time varying 2D field --!
592	IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )'
593	IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
594	IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s'
595	!
596	l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
597	!
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'
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	!
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
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	!
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	!
621	ENDIF
622	!
623	END SUBROUTINE ldf_eiv_init
624
625
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
642	REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars
643	REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace
644	!!----------------------------------------------------------------------
645	!
646	zn (:,:) = 0._wp ! Local initialization
647	zhw(:,:) = 5._wp
648	zah(:,:) = 0._wp
649	zRo(:,:) = 0._wp
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 )
659	zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
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))
663	ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
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 )
677	zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
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))
681	ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
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
688	ENDIF
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 )
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 ) )
695	! Compute aeiw by multiplying Ro^2 and T^-1
696	zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
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.
704	zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
705	zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0
706	END DO
707	END DO
708	CALL lbc_lnk( zaeiw(:,:), 'W', 1. ) ! lateral boundary condition
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
716	CALL lbc_lnk_multi( paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition
717
718	DO jk = 2, jpkm1 !== deeper values equal the surface one ==!
719	paeiu(:,:,jk) = paeiu(:,:,1) * umask(:,:,jk)
720	paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk)
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 ! - -
754	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
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'
761	ENDIF
762
763
764	zpsi_uw(:,:, 1 ) = 0._wp ; zpsi_vw(:,:, 1 ) = 0._wp
765	zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp
766	!
767	DO jk = 2, jpkm1
768	DO jj = 1, jpjm1
769	DO ji = 1, fs_jpim1 ! vector opt.
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)
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
799
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]
812	!
813	INTEGER :: ji, jj, jk ! dummy loop indices
814	REAL(wp) :: zztmp ! local scalar
815	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
816	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
817	!!----------------------------------------------------------------------
818	!
819	!!gm I don't like this routine.... Crazy way of doing things, not optimal at all...
820	!!gm to be redesigned....
821	! !== eiv stream function: output ==!
822	CALL lbc_lnk_multi( psi_uw, 'U', -1. , psi_vw, 'V', -1. )
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	!
829	zw3d(:,:,jpk) = 0._wp ! bottom value always 0
830	!
831	DO jk = 1, jpkm1 ! e2u e3u u_eiv = -dk[psi_uw]
832	zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
833	END DO
834	CALL iom_put( "uoce_eiv", zw3d )
835	!
836	DO jk = 1, jpkm1 ! e1v e3v v_eiv = -dk[psi_vw]
837	zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
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
849	CALL lbc_lnk( zw3d, 'T', 1. ) ! lateral boundary condition
850	CALL iom_put( "woce_eiv", zw3d )
851	!
852	!
853	zztmp = 0.5_wp * rau0 * rcp
854	IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
855	zw2d(:,:) = 0._wp
856	zw3d(:,:,:) = 0._wp
857	DO jk = 1, jpkm1
858	DO jj = 2, jpjm1
859	DO ji = fs_2, fs_jpim1 ! vector opt.
860	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) &
861	& * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) )
862	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
863	END DO
864	END DO
865	END DO
866	CALL lbc_lnk( zw2d, 'U', -1. )
867	CALL lbc_lnk( zw3d, 'U', -1. )
868	CALL iom_put( "ueiv_heattr" , zztmp * zw2d ) ! heat transport in i-direction
869	CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction
870	ENDIF
871	zw2d(:,:) = 0._wp
872	zw3d(:,:,:) = 0._wp
873	DO jk = 1, jpkm1
874	DO jj = 2, jpjm1
875	DO ji = fs_2, fs_jpim1 ! vector opt.
876	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) &
877	& * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) )
878	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
879	END DO
880	END DO
881	END DO
882	CALL lbc_lnk( zw2d, 'V', -1. )
883	CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction
884	CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction
885	!
886	IF( ln_diaptr ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
887	!
888	zztmp = 0.5_wp * 0.5
889	IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
890	zw2d(:,:) = 0._wp
891	zw3d(:,:,:) = 0._wp
892	DO jk = 1, jpkm1
893	DO jj = 2, jpjm1
894	DO ji = fs_2, fs_jpim1 ! vector opt.
895	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) &
896	& * ( tsn (ji,jj,jk,jp_sal) + tsn (ji+1,jj,jk,jp_sal) )
897	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
898	END DO
899	END DO
900	END DO
901	CALL lbc_lnk( zw2d, 'U', -1. )
902	CALL lbc_lnk( zw3d, 'U', -1. )
903	CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction
904	CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction
905	ENDIF
906	zw2d(:,:) = 0._wp
907	zw3d(:,:,:) = 0._wp
908	DO jk = 1, jpkm1
909	DO jj = 2, jpjm1
910	DO ji = fs_2, fs_jpim1 ! vector opt.
911	zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) &
912	& * ( tsn (ji,jj,jk,jp_sal) + tsn (ji,jj+1,jk,jp_sal) )
913	zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
914	END DO
915	END DO
916	END DO
917	CALL lbc_lnk( zw2d, 'V', -1. )
918	CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction
919	CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction
920	!
921	IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
922	!
923	!
924	END SUBROUTINE ldf_eiv_dia
925
926	!!======================================================================
927	END MODULE ldftra

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

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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90 @ 9526

Download in other formats: