Context Navigation

ldfdyn_c2d.h90 @ 7

Last change on this file since 7 was 7, checked in by cholod, 12 years ago
Update of the trunk with the v3.4_r3300
File size: 17.1 KB

Line
1	!!----------------------------------------------------------------------
2	!! * ldfdyn_c2d.h90 *
3	!!----------------------------------------------------------------------
4	!! ldf_dyn_c2d : set the lateral viscosity coefficients
5	!! ldf_dyn_c2d_orca : specific case for orca r2 and r4
6	!!----------------------------------------------------------------------
7
8	!!----------------------------------------------------------------------
9	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
10	!! $Id: ldfdyn_c2d.h90 3294 2012-01-28 16:44:18Z rblod $
11	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
12	!!----------------------------------------------------------------------
13
14	SUBROUTINE ldf_dyn_c2d( ld_print )
15	!!----------------------------------------------------------------------
16	!! * ROUTINE ldf_dyn_c2d *
17	!!
18	!! ** Purpose : initializations of the horizontal ocean physics
19	!!
20	!! ** Method :
21	!! 2D eddy viscosity coefficients ( longitude, latitude )
22	!!
23	!! harmonic operator : ahm1 is defined at t-point
24	!! ahm2 is defined at f-point
25	!! + isopycnal : ahm3 is defined at u-point
26	!! or geopotential ahm4 is defined at v-point
27	!! iso-model level : ahm3, ahm4 not used
28	!!
29	!! biharmonic operator : ahm3 is defined at u-point
30	!! ahm4 is defined at v-point
31	!! : ahm1, ahm2 not used
32	!!
33	!!----------------------------------------------------------------------
34	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
35	!
36	INTEGER :: ji, jj
37	REAL(wp) :: za00, zd_max, zetmax, zeumax, zefmax, zevmax
38	!!----------------------------------------------------------------------
39
40	IF(lwp) WRITE(numout,*)
41	IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
42	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
43
44	! harmonic operator (ahm1, 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	zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
51	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
52
53	IF(lwp) WRITE(numout,*) ' laplacian operator: ahm proportional to e1'
54	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
55
56	za00 = ahm0 / zd_max
57	DO jj = 1, jpj
58	DO ji = 1, jpi
59	zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
60	zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
61	ahm1(ji,jj) = za00 * zetmax
62	ahm2(ji,jj) = za00 * zefmax
63	END DO
64	END DO
65
66	IF( ln_dynldf_iso ) THEN
67	IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum'
68	IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is'
69	IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than'
70	IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.'
71	ENDIF
72
73	! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
74	! ==============================================
75	IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) CALL ldf_dyn_c2d_orca( ld_print )
76	IF( cp_cfg == "orca" .AND. jp_cfg == 1) CALL ldf_dyn_c2d_orca_R1( ld_print )
77
78	! Control print
79	IF( lwp .AND. ld_print ) THEN
80	WRITE(numout,*)
81	WRITE(numout,*) 'inildf: 2D ahm1 array'
82	CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
83	WRITE(numout,*)
84	WRITE(numout,*) 'inildf: 2D ahm2 array'
85	CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
86	ENDIF
87	ENDIF
88
89	! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
90	! ================================= whatever its orientation is)
91	IF( ln_dynldf_bilap ) THEN
92	! (USER: modify ahm3 and ahm4 following your desiderata)
93	! Here: ahm is proportional to the cube of the maximum of the gridspacing
94	! in the to horizontal direction
95
96	zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
97	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
98
99	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
100	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
101
102	za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
103	DO jj = 1, jpj
104	DO ji = 1, jpi
105	zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
106	zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
107	ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
108	ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
109	END DO
110	END DO
111
112	! Control print
113	IF( lwp .AND. ld_print ) THEN
114	WRITE(numout,*)
115	WRITE(numout,*) 'inildf: ahm3 array'
116	CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
117	WRITE(numout,*)
118	WRITE(numout,*) 'inildf: ahm4 array'
119	CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
120	ENDIF
121	ENDIF
122	!
123	END SUBROUTINE ldf_dyn_c2d
124
125
126	SUBROUTINE ldf_dyn_c2d_orca( ld_print )
127	!!----------------------------------------------------------------------
128	!! * ROUTINE ldf_dyn_c2d *
129	!!
130	!! ** W A R N I N G **
131	!!
132	!! ORCA R2 and R4 configurations
133	!!
134	!! ** W A R N I N G **
135	!!
136	!! ** Purpose : initializations of the lateral viscosity for orca R2
137	!!
138	!! ** Method : blah blah blah...
139	!!
140	!!----------------------------------------------------------------------
141	USE ldftra_oce, ONLY: aht0
142	!
143	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
144	!
145	INTEGER :: ji, jj, jn ! dummy loop indices
146	INTEGER :: inum, iim, ijm ! local integers
147	INTEGER :: ifreq, il1, il2, ij, ii
148	REAL(wp) :: zahmeq, zcoft, zcoff, zmsk
149	CHARACTER (len=15) :: clexp
150	INTEGER, POINTER, DIMENSION(:,:) :: icof
151	INTEGER, POINTER, DIMENSION(:,:) :: idata
152	!!----------------------------------------------------------------------
153	!
154	CALL wrk_alloc( jpi , jpj , icof )
155	CALL wrk_alloc( jpidta, jpjdta, idata )
156	!
157	IF(lwp) WRITE(numout,*)
158	IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
159	IF(lwp) WRITE(numout,*) '~~~~~~ --'
160	IF(lwp) WRITE(numout,*) ' orca ocean configuration'
161
162	#if defined key_antarctic
163	# include "ldfdyn_antarctic.h90"
164	#elif defined key_arctic
165	# include "ldfdyn_arctic.h90"
166	#else
167	! Read 2d integer array to specify western boundary increase in the
168	! ===================== equatorial strip (20N-20S) defined at t-points
169
170	CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
171	READ(inum,9101) clexp, iim, ijm
172	READ(inum,'(/)')
173	ifreq = 40
174	il1 = 1
175	DO jn = 1, jpidta/ifreq+1
176	READ(inum,'(/)')
177	il2 = MIN( jpidta, il1+ifreq-1 )
178	READ(inum,9201) ( ii, ji = il1, il2, 5 )
179	READ(inum,'(/)')
180	DO jj = jpjdta, 1, -1
181	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
182	END DO
183	il1 = il1 + ifreq
184	END DO
185
186	DO jj = 1, nlcj
187	DO ji = 1, nlci
188	icof(ji,jj) = idata( mig(ji), mjg(jj) )
189	END DO
190	END DO
191	DO jj = nlcj+1, jpj
192	DO ji = 1, nlci
193	icof(ji,jj) = icof(ji,nlcj)
194	END DO
195	END DO
196	DO jj = 1, jpj
197	DO ji = nlci+1, jpi
198	icof(ji,jj) = icof(nlci,jj)
199	END DO
200	END DO
201
202	9101 FORMAT(1x,a15,2i8)
203	9201 FORMAT(3x,13(i3,12x))
204	9202 FORMAT(i3,41i3)
205
206
207	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operator)
208	! =================
209	! define ahm1 and ahm2 at the right grid point position
210	! (USER: modify ahm1 and ahm2 following your desiderata)
211
212
213	! Decrease ahm to zahmeq m2/s in the tropics
214	! (from 90 to 20 degre: ahm = constant
215	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
216	! from 2.5 to 0 degre: ahm = constant
217	! symmetric in the south hemisphere)
218
219	zahmeq = aht0
220
221	DO jj = 1, jpj
222	DO ji = 1, jpi
223	IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
224	ahm2(ji,jj) = ahm0
225	ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
226	ahm2(ji,jj) = zahmeq
227	ELSE
228	ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2. &
229	* ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
230	ENDIF
231	IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
232	ahm1(ji,jj) = ahm0
233	ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
234	ahm1(ji,jj) = zahmeq
235	ELSE
236	ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2. &
237	* ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
238	ENDIF
239	END DO
240	END DO
241
242	! increase along western boundaries of equatorial strip
243	! t-point
244	DO jj = 1, jpjm1
245	DO ji = 1, jpim1
246	zcoft = FLOAT( icof(ji,jj) ) / 100.
247	ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
248	END DO
249	END DO
250	! f-point
251	icof(:,:) = icof(:,:) * tmask(:,:,1)
252	DO jj = 1, jpjm1
253	DO ji = 1, jpim1 ! NO vector opt.
254	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
255	IF( zmsk == 0. ) THEN
256	zcoff = 1.
257	ELSE
258	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
259	/ (zmsk * 100.)
260	ENDIF
261	ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
262	END DO
263	END DO
264	#endif
265
266	! Lateral boundary conditions on ( ahm1, ahm2 )
267	! ==============
268	CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign
269	CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign
270
271	! Control print
272	IF( lwp .AND. ld_print ) THEN
273	WRITE(numout,*)
274	WRITE(numout,*) 'inildf: 2D ahm1 array'
275	CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
276	WRITE(numout,*)
277	WRITE(numout,*) 'inildf: 2D ahm2 array'
278	CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
279	ENDIF
280	!
281	CALL wrk_dealloc( jpi , jpj , icof )
282	CALL wrk_dealloc( jpidta, jpjdta, idata )
283	!
284	END SUBROUTINE ldf_dyn_c2d_orca
285
286
287	SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
288	!!----------------------------------------------------------------------
289	!! * ROUTINE ldf_dyn_c2d *
290	!!
291	!! ** W A R N I N G **
292	!!
293	!! ORCA R1 configuration
294	!!
295	!! ** W A R N I N G **
296	!!
297	!! ** Purpose : initializations of the lateral viscosity for orca R1
298	!!
299	!! ** Method : blah blah blah...
300	!!
301	!!----------------------------------------------------------------------
302	USE ldftra_oce, ONLY: aht0
303	!
304	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
305	!
306	INTEGER :: ji, jj, jn ! dummy loop indices
307	INTEGER :: inum ! temporary logical unit
308	INTEGER :: iim, ijm
309	INTEGER :: ifreq, il1, il2, ij, ii
310	REAL(wp) :: zahmeq, zcoft, zcoff, zmsk, zam20s
311	CHARACTER (len=15) :: clexp
312	INTEGER, POINTER, DIMENSION(:,:) :: icof
313	INTEGER, POINTER, DIMENSION(:,:) :: idata
314	!!----------------------------------------------------------------------
315	!
316	CALL wrk_alloc( jpi , jpj , icof )
317	CALL wrk_alloc( jpidta, jpjdta, idata )
318	!
319
320	IF(lwp) WRITE(numout,*)
321	IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
322	IF(lwp) WRITE(numout,*) '~~~~~~ --'
323	IF(lwp) WRITE(numout,*) ' orca_r1 configuration'
324
325	#if defined key_antarctic
326	# include "ldfdyn_antarctic.h90"
327	#elif defined key_arctic
328	# include "ldfdyn_arctic.h90"
329	#else
330	! Read 2d integer array to specify western boundary increase in the
331	! ===================== equatorial strip (20N-20S) defined at t-points
332
333	CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', &
334	& 1, numout, lwp )
335	REWIND inum
336	READ(inum,9101) clexp, iim, ijm
337	READ(inum,'(/)')
338	ifreq = 40
339	il1 = 1
340	DO jn = 1, jpidta/ifreq+1
341	READ(inum,'(/)')
342	il2 = MIN( jpidta, il1+ifreq-1 )
343	READ(inum,9201) ( ii, ji = il1, il2, 5 )
344	READ(inum,'(/)')
345	DO jj = jpjdta, 1, -1
346	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
347	END DO
348	il1 = il1 + ifreq
349	END DO
350
351	DO jj = 1, nlcj
352	DO ji = 1, nlci
353	icof(ji,jj) = idata( mig(ji), mjg(jj) )
354	END DO
355	END DO
356	DO jj = nlcj+1, jpj
357	DO ji = 1, nlci
358	icof(ji,jj) = icof(ji,nlcj)
359	END DO
360	END DO
361	DO jj = 1, jpj
362	DO ji = nlci+1, jpi
363	icof(ji,jj) = icof(nlci,jj)
364	END DO
365	END DO
366
367	9101 FORMAT(1x,a15,2i8)
368	9201 FORMAT(3x,13(i3,12x))
369	9202 FORMAT(i3,41i3)
370
371
372	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operator)
373	! =================
374	! define ahm1 and ahm2 at the right grid point position
375	! (USER: modify ahm1 and ahm2 following your desiderata)
376
377
378	! Decrease ahm to zahmeq m2/s in the tropics
379	! (from 90 to 20 degrees: ahm = scaled by local metrics
380	! from 20 to 2.5 degrees: ahm = decrease in (1-cos)/2
381	! from 2.5 to 0 degrees: ahm = constant
382	! symmetric in the south hemisphere)
383
384	zahmeq = aht0
385	zam20s = ahm0COS( rad 20. )
386
387	DO jj = 1, jpj
388	DO ji = 1, jpi
389	IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
390	! leave as set in ldf_dyn_c2d
391	ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
392	ahm2(ji,jj) = zahmeq
393	ELSE
394	ahm2(ji,jj) = zahmeq + (zam20s-zahmeq)/2. &
395	* ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
396	ENDIF
397	IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
398	! leave as set in ldf_dyn_c2d
399	ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
400	ahm1(ji,jj) = zahmeq
401	ELSE
402	ahm1(ji,jj) = zahmeq + (zam20s-zahmeq)/2. &
403	* ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
404	ENDIF
405	END DO
406	END DO
407
408	! increase along western boundaries of equatorial strip
409	! t-point
410	DO jj = 1, jpjm1
411	DO ji = 1, jpim1
412	IF( ABS( gphit(ji,jj) ) < 20. ) THEN
413	zcoft = FLOAT( icof(ji,jj) ) / 100.
414	ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
415	ENDIF
416	END DO
417	END DO
418	! f-point
419	icof(:,:) = icof(:,:) * tmask(:,:,1)
420	DO jj = 1, jpjm1
421	DO ji = 1, jpim1
422	IF( ABS( gphif(ji,jj) ) < 20. ) THEN
423	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
424	IF( zmsk == 0. ) THEN
425	zcoff = 1.
426	ELSE
427	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
428	/ (zmsk * 100.)
429	ENDIF
430	ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
431	ENDIF
432	END DO
433	END DO
434	#endif
435
436	! Lateral boundary conditions on ( ahm1, ahm2 )
437	! ==============
438	CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign
439	CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign
440
441	! Control print
442	IF( lwp .AND. ld_print ) THEN
443	WRITE(numout,*)
444	WRITE(numout,*) 'inildf: 2D ahm1 array'
445	CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
446	WRITE(numout,*)
447	WRITE(numout,*) 'inildf: 2D ahm2 array'
448	CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
449	ENDIF
450	!
451	CALL wrk_dealloc( jpi , jpj , icof )
452	CALL wrk_dealloc( jpidta, jpjdta, idata )
453	!
454	END SUBROUTINE ldf_dyn_c2d_orca_R1

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

Download in other formats:

Original Format