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_c2d.h90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 2633

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

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

Download in other formats:

Original Format