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.

ldfdyn_c3d.h90 in trunk/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 162

Last change on this file since 162 was 162, checked in by opalod, 19 years ago
CT + CL : BUGFIX108 : Add latitudinal variation of viscosity for the ORCA 4� configuration to damp instabilities
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 17.5 KB

Line
1	!!----------------------------------------------------------------------
2	!! * ldfdyn_c3d.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! 'key_dynldf_c3d' 3D lateral eddy viscosity coefficients
7	!!----------------------------------------------------------------------
8
9	SUBROUTINE ldf_dyn_c3d( ld_print )
10	!!----------------------------------------------------------------------
11	!! * ROUTINE ldf_dyn_c3d *
12	!!
13	!! ** Purpose : initializations of the horizontal ocean physics
14	!!
15	!! ** Method : 3D eddy viscosity coef. ( longitude, latitude, depth )
16	!! laplacian operator : ahm1, ahm2 defined at T- and F-points
17	!! ahm2, ahm4 never used
18	!! bilaplacian operator : ahm1, ahm2 never used
19	!! : ahm3, ahm4 defined at U- and V-points
20	!! ??? explanation of the default is missing
21	!!----------------------------------------------------------------------
22	!! * Modules used
23	USE ldftra_oce, ONLY : aht0
24
25	!! * Arguments
26	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
27
28	!! * local variables
29	INTEGER :: ji, jj, jk ! dummy loop indices
30	REAL(wp) :: &
31	zr = 0.2 , & ! maximum of the reduction factor at the bottom ocean
32	! ! ( 0 < zr < 1 )
33	zh = 500., & ! depth of at which start the reduction ( > dept(1) )
34	zdx_max , & ! maximum grid spacing over the global domain
35	za00, zc, zd ! temporary scalars
36	REAL(wp), DIMENSION(jpk) :: zcoef ! temporary workspace
37	!!----------------------------------------------------------------------
38
39	IF(lwp) WRITE(numout,*)
40	IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient'
41	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
42
43
44	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators
45	! ================= whatever its orientation is)
46	IF( ln_dynldf_lap ) THEN
47	! define ahm1 and ahm2 at the right grid point position
48	! (USER: modify ahm1 and ahm2 following your desiderata)
49
50	zdx_max = MAXVAL( e1t(:,:) )
51	IF( lk_mpp ) CALL mpp_max( zdx_max ) ! max over the global domain
52
53	IF(lwp) WRITE(numout,*) ' laplacian operator: ahm proportional to e1'
54	IF(lwp) WRITE(numout,*) ' Caution, here we assume your mesh is isotropic ...'
55	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
56
57
58	za00 = ahm0 / zdx_max
59
60	IF( ln_dynldf_iso ) THEN
61	IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum'
62	IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is'
63	IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than'
64	IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.'
65	ENDIF
66
67	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 ) ! vertical profile
68	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 ) ! vertical profile
69	DO jk = 1,jpk
70	ahm1(:,:,jk) = za00 * e1t(:,:) * ahm1(:,:,jk)
71	ahm2(:,:,jk) = za00 * e1f(:,:) * ahm2(:,:,jk)
72	END DO
73
74
75	! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2)
76	! ==============================================
77	IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
78	IF(lwp) WRITE(numout,*)
79	IF(lwp) WRITE(numout,*) ' ORCA R2 or R4: overwrite the previous definition of ahm'
80	IF(lwp) WRITE(numout,*) ' ============='
81	CALL ldf_dyn_c3d_orca( ld_print )
82	ENDIF
83
84	ENDIF
85
86	! Control print
87	IF(lwp .AND. ld_print ) THEN
88	WRITE(numout,*)
89	WRITE(numout,*) ' 3D ahm1 array (k=1)'
90	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
91	WRITE(numout,*)
92	WRITE(numout,*) ' 3D ahm2 array (k=1)'
93	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
94	ENDIF
95
96
97	! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
98	! ================================ whatever its orientation is)
99	! (USER: modify ahm3 and ahm4 following your desiderata)
100	! Here: ahm is proportional to the cube of the maximum of the gridspacing
101	! in the to horizontal direction
102
103	IF( ln_dynldf_bilap ) THEN
104
105	zdx_max = MAXVAL( e1u(:,:) )
106	IF( lk_mpp ) CALL mpp_max( zdx_max ) ! max over the global domain
107
108	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
109	IF(lwp) WRITE(numout,*) ' Caution, here we assume your mesh is isotropic ...'
110	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
111
112	za00 = ahm0 / ( zdx_max * zdx_max * zdx_max )
113	ahm3(:,:,1) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:)
114	ahm4(:,:,1) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:)
115
116	zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0
117	IF( lk_zco ) THEN ! z-coordinate, same profile everywhere
118	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
119	DO jk = 1, jpk
120	IF( fsdept(1,1,jk) <= zh ) THEN
121	zcoef(jk) = 1.e0
122	ELSE
123	zcoef(jk) = 1.e0 + ( zr - 1.e0 ) &
124	& * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) &
125	& / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) )
126	ENDIF
127	ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
128	ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
129	IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
130	END DO
131	ELSE ! partial steps or s-ccordinate
132	# if defined key_partial_steps \|\| defined key_s_coord
133	zc = MAXVAL( fsdept(:,:,jpkm1) )
134	# else
135	zc = fsdept(:,:,jpkm1)
136	# endif
137	IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain
138
139	zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) )
140	DO jk = 2, jpkm1
141	DO jj = 1, jpj
142	DO ji = 1, jpi
143	IF( fsdept(ji,jj,jk) <= zh ) THEN
144	ahm3(ji,jj,jk) = ahm3(ji,jj,1)
145	ahm4(ji,jj,jk) = ahm4(ji,jj,1)
146	ELSE
147	zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc
148	ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
149	ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
150	ENDIF
151	END DO
152	END DO
153	END DO
154	ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
155	ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
156	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
157	DO jk = 1, jpk
158	IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
159	END DO
160	ENDIF
161
162	! Control print
163	IF( lwp .AND. ld_print ) THEN
164	WRITE(numout,*)
165	WRITE(numout,*) 'inildf: ahm3 array at level 1'
166	CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
167	WRITE(numout,*)
168	WRITE(numout,*) 'inildf: ahm4 array at level 1'
169	CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
170	ENDIF
171	ENDIF
172
173	END SUBROUTINE ldf_dyn_c3d
174
175
176	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
177	!!----------------------------------------------------------------------
178	!! * ROUTINE ldf_dyn_c3d *
179	!!
180	!! ** Purpose : ORCA R2 an R4 only
181	!!
182	!! ** Method : blah blah blah ....
183	!!----------------------------------------------------------------------
184	!! * Modules used
185	USE ldftra_oce, ONLY : aht0
186
187	!! * Arguments
188	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
189
190	!! * local variables
191	INTEGER :: ji, jj, jk, jn ! dummy loop indices
192	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
193	INTEGER :: inum = 11 ! temporary logical unit
194	INTEGER :: iost, iim, ijm
195	INTEGER :: ifreq, il1, il2, ij, ii
196	INTEGER, DIMENSION(jpidta, jpjdta) :: idata
197	INTEGER, DIMENSION(jpi , jpj ) :: icof
198
199	REAL(wp) :: &
200	zahmeq, zcoff, zcoft, zmsk, & ! ???
201	zemax, zemin, zeref, zahmm
202	REAL(wp), DIMENSION(jpi,jpj) :: zahm0
203	REAL(wp), DIMENSION(jpk) :: zcoef
204
205	CHARACTER (len=15) :: clexp
206	!!----------------------------------------------------------------------
207
208	IF(lwp) WRITE(numout,*)
209	IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
210	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
211	IF(lwp) WRITE(numout,*)
212	IF(lwp) WRITE(numout,*) ' orca R2 or R4 ocean model'
213	IF(lwp) WRITE(numout,*) ' reduced in the surface Eq. strip '
214	IF(lwp) WRITE(numout,*)
215
216	! Read 2d integer array to specify western boundary increase in the
217	! ===================== equatorial strip (20N-20S) defined at t-points
218
219	OPEN( UNIT=inum, FILE='ahmcoef', STATUS='OLD', &
220	& FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 , &
221	& IOSTAT= iost)
222	IF( iost == 0 ) THEN
223	IF(lwp) THEN
224	WRITE(numout,*) ' file : ahmcoef open ok'
225	WRITE(numout,*) ' unit = ', inum
226	WRITE(numout,*) ' status = OLD'
227	WRITE(numout,*) ' form = FORMATTED'
228	WRITE(numout,*) ' access = SEQUENTIAL'
229	WRITE(numout,*)
230	ENDIF
231	ENDIF
232	111 CONTINUE
233	IF( iost /= 0 ) THEN
234	IF(lwp) THEN
235	WRITE(numout,*)
236	WRITE(numout,*) ' ===>>>> : bad opening file: ahmcoef, we stop. verify the file '
237	WRITE(numout,*) ' ======= === '
238	ENDIF
239	nstop = nstop + 1
240	ENDIF
241
242	REWIND inum
243	READ(inum,9101) clexp, iim, ijm
244	READ(inum,'(/)')
245	ifreq = 40
246	il1 = 1
247	DO jn = 1, jpidta/ifreq+1
248	READ(inum,'(/)')
249	il2 = MIN( jpidta, il1+ifreq-1 )
250	READ(inum,9201) ( ii, ji = il1, il2, 5 )
251	READ(inum,'(/)')
252	DO jj = jpjdta, 1, -1
253	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
254	END DO
255	il1 = il1 + ifreq
256	END DO
257
258	DO jj = 1, nlcj
259	DO ji = 1, nlci
260	icof(ji,jj) = idata( mig(ji), mjg(jj) )
261	END DO
262	END DO
263	DO jj = nlcj+1, jpj
264	DO ji = 1, nlci
265	icof(ji,jj) = icof(ji,nlcj)
266	END DO
267	END DO
268	DO jj = 1, jpj
269	DO ji = nlci+1, jpi
270	icof(ji,jj) = icof(nlci,jj)
271	END DO
272	END DO
273
274	9101 FORMAT(1x,a15,2i8)
275	9201 FORMAT(3x,13(i3,12x))
276	9202 FORMAT(i3,41i3)
277
278	! Set ahm1 and ahm2
279	! =================
280
281	! define ahm1 and ahm2 at the right grid point position
282	! (USER: modify ahm1 and ahm2 following your desiderata)
283	! biharmonic : ahm1 (ahm2) defined at u- (v-) point
284	! harmonic : ahm1 (ahm2) defined at t- (f-) point
285
286	! first level : as for 2D coefficients
287
288	! Decrease ahm to zahmeq m2/s in the tropics
289	! (from 90 to 20 degre: ahm = constant
290	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
291	! from 2.5 to 0 degre: ahm = constant
292	! symmetric in the south hemisphere)
293
294	IF( jp_cfg == 4 ) THEN
295	zahmeq = 5.0 * aht0
296	zahmm = min( 160000.0, ahm0)
297	zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
298	zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
299	zeref = MAXVAL ( e1t(:,:) * e2t(:,:), &
300	& tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
301
302	DO jj = 1, jpj
303	DO ji = 1, jpi
304	zmsk = e1t(ji,jj) * e2t(ji,jj)
305	IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
306	zahm0(ji,jj) = ahm0
307	ELSE
308	IF( zmsk .GE. zeref ) THEN
309	zahm0(ji,jj) = ahm0
310	ELSE
311	zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - &
312	& cos((rpi0.5(zmsk-zemin)/(zeref-zemin))))
313	ENDIF
314	ENDIF
315	END DO
316	END DO
317	ENDIF
318
319	IF( jp_cfg == 2 ) THEN
320	zahmeq = aht0
321	zahmm = ahm0
322	zahm0(:,:) = ahm0
323	ENDIF
324
325	DO jj = 1, jpj
326	DO ji = 1, jpi
327	IF( ABS(gphif(ji,jj)) >= 20.) THEN
328	ahm2(ji,jj,1) = zahm0(ji,jj)
329	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
330	ahm2(ji,jj,1) = zahmeq
331	ELSE
332	ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
333	& (1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
334	ENDIF
335	IF( ABS(gphit(ji,jj)) >= 20.) THEN
336	ahm1(ji,jj,1) = zahm0(ji,jj)
337	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
338	ahm1(ji,jj,1) = zahmeq
339	ELSE
340	ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
341	& (1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
342	ENDIF
343	END DO
344	END DO
345
346	! increase along western boundaries of equatorial strip
347	! t-point
348	DO jj = 1, jpjm1
349	DO ji = 1, jpim1
350	zcoft = float( icof(ji,jj) ) / 100.
351	ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
352	END DO
353	END DO
354	! f-point
355	icof(:,:) = icof(:,:) * tmask(:,:,1)
356	DO jj = 1, jpjm1
357	DO ji = 1, jpim1
358	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
359	IF( zmsk == 0. ) THEN
360	zcoff = 1.
361	ELSE
362	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
363	/ (zmsk * 100.)
364	ENDIF
365	ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
366	END DO
367	END DO
368
369	! other level: re-increase the coef in the deep ocean
370
371	DO jk = 1, 21
372	zcoef(jk) = 1.
373	END DO
374	zcoef(22) = 2.
375	zcoef(23) = 3.
376	zcoef(24) = 5.
377	zcoef(25) = 7.
378	zcoef(26) = 9.
379	DO jk = 27, jpk
380	zcoef(jk) = 10.
381	END DO
382
383	DO jk = 2, jpk
384	ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
385	ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
386	END DO
387
388	IF( jp_cfg == 4 ) THEN
389	! Limit AHM in Gibraltar strait
390	ij0 = 50 ; ij1 = 53
391	ii0 = 69 ; ii1 = 71
392	DO jk = 1, jpk
393	ahm1(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm1 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) )
394	ahm2(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm2 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) )
395	END DO
396	ENDIF
397
398	! Lateral boundary conditions on ( ahm1, ahm2 )
399	! ==============
400	CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign
401	CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign
402
403	! Control print
404
405	IF(lwp) THEN
406	WRITE(numout,*)
407	WRITE(numout,*) ' 3D ahm1 array (k=1)'
408	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
409	WRITE(numout,*)
410	WRITE(numout,*) ' 3D ahm2 array (k=1)'
411	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
412	WRITE(numout,*)
413	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
414	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
415	ENDIF
416
417
418	! Set ahm3 and ahm4
419	! =================
420
421	! define ahm3 and ahm4 at the right grid point position
422	! initialization to a constant value
423	! (USER: modify ahm3 and ahm4 following your desiderata)
424	! harmonic isopycnal or geopotential:
425	! ahm3 (ahm4) defined at u- (v-) point
426	DO jk = 1, jpk
427	DO jj = 2, jpj
428	DO ji = 2, jpi
429	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
430	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
431	END DO
432	END DO
433	END DO
434	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
435	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
436
437	! Lateral boundary conditions on ( ahm3, ahm4 )
438	! ==============
439	CALL lbc_lnk( ahm3, 'U', 1. ) ! U-point, unchanged sign
440	CALL lbc_lnk( ahm4, 'V', 1. ) ! V-point, unchanged sign
441
442	! Control print
443
444	IF( lwp .AND. ld_print ) THEN
445	WRITE(numout,*)
446	WRITE(numout,*) ' ahm3 array level 1'
447	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
448	WRITE(numout,*)
449	WRITE(numout,*) ' ahm4 array level 1'
450	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
451	ENDIF
452
453	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format