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

Context Navigation

source: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 247

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

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

Download in other formats:

Original Format