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 tags/nemo_dev_x3/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

source: tags/nemo_dev_x3/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 105

Last change on this file since 105 was 105, checked in by cvs2svn, 20 years ago
This commit was manufactured by cvs2svn to create tag 'nemo_dev_x3'.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 11.4 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	!! OPA 9.0 , LODYC-IPSL (2003)
10	!!----------------------------------------------------------------------
11
12	SUBROUTINE ldf_dyn_c2d( ld_print )
13	!!----------------------------------------------------------------------
14	!! * ROUTINE ldf_dyn_c2d *
15	!!
16	!! ** Purpose : initializations of the horizontal ocean physics
17	!!
18	!! ** Method :
19	!! 2D eddy viscosity coefficients ( longitude, latitude )
20	!!
21	!! harmonic operator : ahm1 is defined at t-point
22	!! ahm2 is defined at f-point
23	!! + isopycnal : ahm3 is defined at u-point
24	!! or geopotential ahm4 is defined at v-point
25	!! iso-model level : ahm3, ahm4 not used
26	!!
27	!! biharmonic operator : ahm1 is defined at u-point
28	!! ahm2 is defined at v-point
29	!! : ahm3, ahm4 not used
30	!!
31	!!----------------------------------------------------------------------
32	!! * Arguments
33	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
34
35	!! * Local variables
36	REAL(wp) :: za00, zdx_max
37	!!----------------------------------------------------------------------
38
39	IF(lwp) WRITE(numout,*)
40	IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
41	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
42	IF(lwp) WRITE(numout,*)
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	za00 = ahm0 / zdx_max
58	ahm1(:,:) = za00 * e1t(:,:)
59	ahm2(:,:) = za00 * e1f(:,:)
60
61	IF( ln_dynldf_iso ) THEN
62	IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum'
63	IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is'
64	IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than'
65	IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.'
66	ENDIF
67
68	! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2)
69	! ==============================================
70	IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) CALL ldf_dyn_c2d_orca( ld_print )
71
72	! Control print
73	IF( lwp .AND. ld_print ) THEN
74	WRITE(numout,*)
75	WRITE(numout,*) 'inildf: 2D ahm1 array'
76	CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
77	WRITE(numout,*)
78	WRITE(numout,*) 'inildf: 2D ahm2 array'
79	CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
80	ENDIF
81	ENDIF
82
83	! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
84	! ================================ whatever its orientation is)
85	IF( ln_dynldf_bilap ) THEN
86	! (USER: modify ahm3 and ahm4 following your desiderata)
87	! Here: ahm is proportional to the cube of the maximum of the gridspacing
88	! in the to horizontal direction
89
90	zdx_max = MAXVAL( e1u(:,:) )
91	IF( lk_mpp ) CALL mpp_max( zdx_max ) ! max over the global domain
92
93	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
94	IF(lwp) WRITE(numout,*) ' Caution, here we assume your mesh is isotropic ...'
95	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
96
97	za00 = ahm0 / ( zdx_max * zdx_max * zdx_max )
98	ahm3(:,:) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:)
99	ahm4(:,:) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:)
100
101	! Control print
102	IF( lwp .AND. ld_print ) THEN
103	WRITE(numout,*)
104	WRITE(numout,*) 'inildf: ahm3 array'
105	CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
106	WRITE(numout,*)
107	WRITE(numout,*) 'inildf: ahm4 array'
108	CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
109	ENDIF
110	ENDIF
111
112
113	END SUBROUTINE ldf_dyn_c2d
114
115
116	SUBROUTINE ldf_dyn_c2d_orca( ld_print )
117	!!----------------------------------------------------------------------
118	!! * ROUTINE ldf_dyn_c2d *
119	!!
120	!! ** W A R N I N G **
121	!!
122	!! ORCA R2 and R4 configurations
123	!!
124	!! ** W A R N I N G **
125	!!
126	!! ** Purpose : initializations of the lateral viscosity for orca R2
127	!!
128	!! ** Method : blah blah blah...
129	!!
130	!!----------------------------------------------------------------------
131	!! * Modules used
132	USE ldftra_oce, ONLY : aht0
133
134	!! * Arguments
135	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
136
137	!! * Local variables
138	INTEGER :: ji, jj, jn ! dummy loop indices
139	INTEGER :: inum = 11 ! temporary logical unit
140	INTEGER :: iost, iim, ijm
141	INTEGER :: ifreq, il1, il2, ij, ii
142	INTEGER, DIMENSION(jpidta,jpidta) :: idata
143	INTEGER, DIMENSION(jpi ,jpj ) :: icof
144
145	REAL(wp) :: zahmeq, zcoft, zcoff, zmsk
146
147	CHARACTER (len=15) :: clexp
148	!!----------------------------------------------------------------------
149
150	IF(lwp) WRITE(numout,*)
151	IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
152	IF(lwp) WRITE(numout,*) '~~~~~~ --'
153	IF(lwp) WRITE(numout,*)
154	IF(lwp) WRITE(numout,*) ' orca ocean model'
155	IF(lwp) WRITE(numout,*)
156
157	#if defined key_antarctic
158	# include "ldfdyn_antarctic.h90"
159	#elif defined key_arctic
160	# include "ldfdyn_arctic.h90"
161	#else
162	! Read 2d integer array to specify western boundary increase in the
163	! ===================== equatorial strip (20N-20S) defined at t-points
164
165	OPEN( UNIT=inum, FILE='ahmcoef', STATUS='OLD', &
166	& FORM='FORMATTED', ACCESS='SEQUENTIAL', ERR=111 , &
167	& IOSTAT= iost )
168	IF( iost == 0 ) THEN
169	IF(lwp) WRITE(numout,*) ' file : ahmcoef open ok'
170	IF(lwp) WRITE(numout,*) ' unit = ', inum
171	IF(lwp) WRITE(numout,*) ' status = OLD'
172	IF(lwp) WRITE(numout,*) ' form = FORMATTED'
173	IF(lwp) WRITE(numout,*) ' access = SEQUENTIAL'
174	IF(lwp) WRITE(numout,*)
175	ENDIF
176	111 CONTINUE
177	IF( iost /= 0 ) THEN
178	IF(lwp) WRITE(numout,*)
179	IF(lwp) WRITE(numout,*) ' ===>>>> : bad opening file: ahmcoef'
180	IF(lwp) WRITE(numout,*) ' ======= === '
181	IF(lwp) WRITE(numout,*) ' we stop. verify the file '
182	IF(lwp) WRITE(numout,*)
183	STOP 'ldfdyn_c2d.h90'
184	ENDIF
185
186	REWIND inum
187	READ(inum,9101) clexp, iim, ijm
188	READ(inum,'(/)')
189	ifreq = 40
190	il1 = 1
191	DO jn = 1, jpidta/ifreq+1
192	READ(inum,'(/)')
193	il2 = MIN( jpidta, il1+ifreq-1 )
194	READ(inum,9201) ( ii, ji = il1, il2, 5 )
195	READ(inum,'(/)')
196	DO jj = jpjdta, 1, -1
197	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
198	END DO
199	il1 = il1 + ifreq
200	END DO
201
202	DO jj = 1, nlcj
203	DO ji = 1, nlci
204	icof(ji,jj) = idata( mig(ji), mjg(jj) )
205	END DO
206	END DO
207	DO jj = nlcj+1, jpj
208	DO ji = 1, nlci
209	icof(ji,jj) = icof(ji,nlcj)
210	END DO
211	END DO
212	DO jj = 1, jpj
213	DO ji = nlci+1, jpi
214	icof(ji,jj) = icof(nlci,jj)
215	END DO
216	END DO
217
218	9101 FORMAT(1x,a15,2i8)
219	9201 FORMAT(3x,13(i3,12x))
220	9202 FORMAT(i3,41i3)
221
222
223	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operator)
224	! =================
225	! define ahm1 and ahm2 at the right grid point position
226	! (USER: modify ahm1 and ahm2 following your desiderata)
227
228
229	! Decrease ahm to zahmeq m2/s in the tropics
230	! (from 90 to 20 degre: ahm = constant
231	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
232	! from 2.5 to 0 degre: ahm = constant
233	! symmetric in the south hemisphere)
234
235	zahmeq = aht0
236
237	DO jj = 1, jpj
238	DO ji = 1, jpi
239	IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
240	ahm2(ji,jj) = ahm0
241	ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
242	ahm2(ji,jj) = zahmeq
243	ELSE
244	ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2. &
245	* ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
246	ENDIF
247	IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
248	ahm1(ji,jj) = ahm0
249	ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
250	ahm1(ji,jj) = zahmeq
251	ELSE
252	ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2. &
253	* ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
254	ENDIF
255	END DO
256	END DO
257
258	! increase along western boundaries of equatorial strip
259	! t-point
260	DO jj = 1, jpjm1
261	DO ji = 1, jpim1
262	zcoft = FLOAT( icof(ji,jj) ) / 100.
263	ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
264	END DO
265	END DO
266	! f-point
267	icof(:,:) = icof(:,:) * tmask(:,:,1)
268	DO jj = 1, jpjm1
269	DO ji = 1, jpim1
270	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
271	IF( zmsk == 0. ) THEN
272	zcoff = 1.
273	ELSE
274	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
275	/ (zmsk * 100.)
276	ENDIF
277	ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
278	END DO
279	END DO
280	#endif
281
282	! Lateral boundary conditions on ( ahm1, ahm2 )
283	! ==============
284	CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign
285	CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign
286
287	! Control print
288	IF( lwp .AND. ld_print ) THEN
289	WRITE(numout,*)
290	WRITE(numout,*) 'inildf: 2D ahm1 array'
291	CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
292	WRITE(numout,*)
293	WRITE(numout,*) 'inildf: 2D ahm2 array'
294	CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
295	ENDIF
296
297	END SUBROUTINE ldf_dyn_c2d_orca

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

Download in other formats:

Original Format