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

Last change on this file since 32 was 32, checked in by opalod, 20 years ago
CT : UPDATE001 : First major NEMO update
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 15.8 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	zc = MAXVAL( fsdept(:,:,jpkm1) )
133	IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain
134
135	zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) )
136	DO jk = 2, jpkm1
137	DO jj = 1, jpj
138	DO ji = 1, jpi
139	IF( fsdept(ji,jj,jk) <= zh ) THEN
140	ahm3(ji,jj,jk) = ahm3(ji,jj,1)
141	ahm4(ji,jj,jk) = ahm4(ji,jj,1)
142	ELSE
143	zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc
144	ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
145	ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
146	ENDIF
147	END DO
148	END DO
149	END DO
150	ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
151	ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
152	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
153	DO jk = 1, jpk
154	IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
155	END DO
156	ENDIF
157
158	! Control print
159	IF( lwp .AND. ld_print ) THEN
160	WRITE(numout,*)
161	WRITE(numout,*) 'inildf: ahm3 array at level 1'
162	CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
163	WRITE(numout,*)
164	WRITE(numout,*) 'inildf: ahm4 array at level 1'
165	CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
166	ENDIF
167	ENDIF
168
169	END SUBROUTINE ldf_dyn_c3d
170
171
172	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
173	!!----------------------------------------------------------------------
174	!! * ROUTINE ldf_dyn_c3d *
175	!!
176	!! ** Purpose : ORCA R2 an R4 only
177	!!
178	!! ** Method : blah blah blah ....
179	!!----------------------------------------------------------------------
180	!! * Modules used
181	USE ldftra_oce, ONLY : aht0
182
183	!! * Arguments
184	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
185
186	!! * local variables
187	INTEGER :: ji, jj, jk, jn ! dummy loop indices
188	INTEGER :: inum = 11 ! temporary logical unit
189	INTEGER :: iost, iim, ijm
190	INTEGER :: ifreq, il1, il2, ij, ii
191	INTEGER, DIMENSION(jpidta, jpjdta) :: idata
192	INTEGER, DIMENSION(jpi , jpj ) :: icof
193
194	REAL(wp) :: zahmeq, zcoff, zcoft, zmsk
195	REAL(wp), DIMENSION(jpk) :: zcoef
196
197	CHARACTER (len=15) :: clexp
198	!!----------------------------------------------------------------------
199
200	IF(lwp) WRITE(numout,*)
201	IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
202	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
203	IF(lwp) WRITE(numout,*)
204	IF(lwp) WRITE(numout,*) ' orca R2 or R4 ocean model'
205	IF(lwp) WRITE(numout,*) ' reduced in the surface Eq. strip '
206	IF(lwp) WRITE(numout,*)
207
208	! Read 2d integer array to specify western boundary increase in the
209	! ===================== equatorial strip (20N-20S) defined at t-points
210
211	OPEN( UNIT=inum, FILE='ahmcoef', STATUS='OLD', &
212	& FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 , &
213	& IOSTAT= iost)
214	IF( iost == 0 ) THEN
215	IF(lwp) THEN
216	WRITE(numout,*) ' file : ahmcoef open ok'
217	WRITE(numout,*) ' unit = ', inum
218	WRITE(numout,*) ' status = OLD'
219	WRITE(numout,*) ' form = FORMATTED'
220	WRITE(numout,*) ' access = SEQUENTIAL'
221	WRITE(numout,*)
222	ENDIF
223	ENDIF
224	111 CONTINUE
225	IF( iost /= 0 ) THEN
226	IF(lwp) THEN
227	WRITE(numout,*)
228	WRITE(numout,*) ' ===>>>> : bad opening file: ahmcoef, we stop. verify the file '
229	WRITE(numout,*) ' ======= === '
230	ENDIF
231	nstop = nstop + 1
232	ENDIF
233
234	REWIND inum
235	READ(inum,9101) clexp, iim, ijm
236	READ(inum,'(/)')
237	ifreq = 40
238	il1 = 1
239	DO jn = 1, jpidta/ifreq+1
240	READ(inum,'(/)')
241	il2 = MIN( jpidta, il1+ifreq-1 )
242	READ(inum,9201) ( ii, ji = il1, il2, 5 )
243	READ(inum,'(/)')
244	DO jj = jpjdta, 1, -1
245	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
246	END DO
247	il1 = il1 + ifreq
248	END DO
249
250	DO jj = 1, nlcj
251	DO ji = 1, nlci
252	icof(ji,jj) = idata( mig(ji), mjg(jj) )
253	END DO
254	END DO
255	DO jj = nlcj+1, jpj
256	DO ji = 1, nlci
257	icof(ji,jj) = icof(ji,nlcj)
258	END DO
259	END DO
260	DO jj = 1, jpj
261	DO ji = nlci+1, jpi
262	icof(ji,jj) = icof(nlci,jj)
263	END DO
264	END DO
265
266	9101 FORMAT(1x,a15,2i8)
267	9201 FORMAT(3x,13(i3,12x))
268	9202 FORMAT(i3,41i3)
269
270	! Set ahm1 and ahm2
271	! =================
272
273	! define ahm1 and ahm2 at the right grid point position
274	! (USER: modify ahm1 and ahm2 following your desiderata)
275	! biharmonic : ahm1 (ahm2) defined at u- (v-) point
276	! harmonic : ahm1 (ahm2) defined at t- (f-) point
277
278	! first level : as for 2D coefficients
279
280	! Decrease ahm to zahmeq m2/s in the tropics
281	! (from 90 to 20 degre: ahm = constant
282	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
283	! from 2.5 to 0 degre: ahm = constant
284	! symmetric in the south hemisphere)
285
286	IF( jp_cfg == 4 ) zahmeq = 5.0 * aht0
287	IF( jp_cfg == 2 ) zahmeq = aht0
288
289	DO jj = 1, jpj
290	DO ji = 1, jpi
291	IF( ABS(gphif(ji,jj)) >= 20.) THEN
292	ahm2(ji,jj,1) = ahm0
293	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
294	ahm2(ji,jj,1) = zahmeq
295	ELSE
296	ahm2(ji,jj,1) = zahmeq &
297	+ (ahm0-zahmeq)/2.(1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
298	ENDIF
299	IF( ABS(gphit(ji,jj)) >= 20.) THEN
300	ahm1(ji,jj,1) = ahm0
301	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
302	ahm1(ji,jj,1) = zahmeq
303	ELSE
304	ahm1(ji,jj,1) = zahmeq &
305	+ (ahm0-zahmeq)/2.(1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
306	ENDIF
307	END DO
308	END DO
309
310	! increase along western boundaries of equatorial strip
311	! t-point
312	DO jj = 1, jpjm1
313	DO ji = 1, jpim1
314	zcoft = float( icof(ji,jj) ) / 100.
315	ahm1(ji,jj,1) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj,1)
316	END DO
317	END DO
318	! f-point
319	icof(:,:) = icof(:,:) * tmask(:,:,1)
320	DO jj = 1, jpjm1
321	DO ji = 1, jpim1
322	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
323	IF( zmsk == 0. ) THEN
324	zcoff = 1.
325	ELSE
326	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
327	/ (zmsk * 100.)
328	ENDIF
329	ahm2(ji,jj,1) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj,1)
330	END DO
331	END DO
332
333	! other level: re-increase the coef in the deep ocean
334
335	DO jk = 1, 21
336	zcoef(jk) = 1.
337	END DO
338	zcoef(22) = 2.
339	zcoef(23) = 3.
340	zcoef(24) = 5.
341	zcoef(25) = 7.
342	zcoef(26) = 9.
343	DO jk = 27, jpk
344	zcoef(jk) = 10.
345	END DO
346
347	DO jk = 2, jpk
348	ahm1(:,:,jk) = MIN( ahm0, zcoef(jk) * ahm1(:,:,1) )
349	ahm2(:,:,jk) = MIN( ahm0, zcoef(jk) * ahm2(:,:,1) )
350	END DO
351
352	! Lateral boundary conditions on ( ahm1, ahm2 )
353	! ==============
354	CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign
355	CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign
356
357	! Control print
358
359	IF(lwp) THEN
360	WRITE(numout,*)
361	WRITE(numout,*) ' 3D ahm1 array (k=1)'
362	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
363	WRITE(numout,*)
364	WRITE(numout,*) ' 3D ahm2 array (k=1)'
365	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
366	WRITE(numout,*)
367	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
368	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
369	ENDIF
370
371
372	! Set ahm3 and ahm4
373	! =================
374
375	! define ahm3 and ahm4 at the right grid point position
376	! initialization to a constant value
377	! (USER: modify ahm3 and ahm4 following your desiderata)
378	! harmonic isopycnal or geopotential:
379	! ahm3 (ahm4) defined at u- (v-) point
380	DO jk = 1, jpk
381	DO jj = 2, jpj
382	DO ji = 2, jpi
383	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
384	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
385	END DO
386	END DO
387	END DO
388	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
389	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
390
391	! Lateral boundary conditions on ( ahm3, ahm4 )
392	! ==============
393	CALL lbc_lnk( ahm3, 'U', 1. ) ! U-point, unchanged sign
394	CALL lbc_lnk( ahm4, 'V', 1. ) ! V-point, unchanged sign
395
396	! Control print
397
398	IF( lwp .AND. ld_print ) THEN
399	WRITE(numout,*)
400	WRITE(numout,*) ' ahm3 array level 1'
401	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
402	WRITE(numout,*)
403	WRITE(numout,*) ' ahm4 array level 1'
404	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
405	ENDIF
406
407	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format