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 @ 901

Last change on this file since 901 was 901, checked in by rblod, 16 years ago
Change ldf computation to take into account the grid anisotropy, see ticket #115
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 18.0 KB

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

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

Download in other formats:

Original Format