Context Navigation

ldfdyn_c3d.h90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 18.3 KB

Line
1	!!----------------------------------------------------------------------
2	!! * ldfdyn_c3d.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
7	!! $Id$
8	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
9	!!----------------------------------------------------------------------
10
11	!!----------------------------------------------------------------------
12	!! 'key_dynldf_c3d' 3D lateral eddy viscosity coefficients
13	!!----------------------------------------------------------------------
14
15	SUBROUTINE ldf_dyn_c3d( ld_print )
16	!!----------------------------------------------------------------------
17	!! * ROUTINE ldf_dyn_c3d *
18	!!
19	!! ** Purpose : initializations of the horizontal ocean physics
20	!!
21	!! ** Method : 3D eddy viscosity coef. ( longitude, latitude, depth )
22	!! laplacian operator : ahm1, ahm2 defined at T- and F-points
23	!! ahm2, ahm4 never used
24	!! bilaplacian operator : ahm1, ahm2 never used
25	!! : ahm3, ahm4 defined at U- and V-points
26	!! ??? explanation of the default is missing
27	!!----------------------------------------------------------------------
28	USE ldftra_oce, ONLY : aht0
29	USE iom
30	!!
31	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
32	!!
33	INTEGER :: ji, jj, jk ! dummy loop indices
34	REAL(wp) :: zr = 0.2 ! maximum of the reduction factor at the bottom ocean ( 0 < zr < 1 )
35	REAL(wp) :: zh = 500. ! depth of at which start the reduction ( > dept(1) )
36	REAL(wp) :: zd_max ! maximum grid spacing over the global domain
37	REAL(wp) :: za00, zc, zd, zetmax, zefmax, zeumax, zevmax ! local scalars
38	REAL(wp), POINTER, DIMENSION(:) :: zcoef
39	!!----------------------------------------------------------------------
40	!
41	CALL wrk_alloc( jpk, zcoef )
42	!
43	IF(lwp) WRITE(numout,*)
44	IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient'
45	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
46
47
48	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators
49	! ================= whatever its orientation is)
50	IF( ln_dynldf_lap ) THEN
51	! define ahm1 and ahm2 at the right grid point position
52	! (USER: modify ahm1 and ahm2 following your desiderata)
53
54	zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
55	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
56
57	IF(lwp) WRITE(numout,*) ' laplacian operator: ahm proportional to e1'
58	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
59
60	za00 = ahm0 / zd_max
61
62	IF( ln_dynldf_iso ) THEN
63	IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum'
64	IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is'
65	IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than'
66	IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.'
67	ENDIF
68
69	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 ) ! vertical profile
70	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 ) ! vertical profile
71	DO jk = 1,jpk
72	DO jj = 1, jpj
73	DO ji = 1, jpi
74	zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
75	zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
76	ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk)
77	ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk)
78	END DO
79	END DO
80	END DO
81
82
83	! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
84	! ==============================================
85	IF( cp_cfg == "orca" .AND. ( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
86	IF(lwp) WRITE(numout,*)
87	IF(lwp) WRITE(numout,*) ' ORCA R1, R2 or R4: overwrite the previous definition of ahm'
88	IF(lwp) WRITE(numout,*) ' ================='
89	CALL ldf_dyn_c3d_orca( ld_print )
90	ENDIF
91
92	ENDIF
93
94	! Control print
95	IF(lwp .AND. ld_print ) THEN
96	WRITE(numout,*)
97	WRITE(numout,*) ' 3D ahm1 array (k=1)'
98	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
99	WRITE(numout,*)
100	WRITE(numout,*) ' 3D ahm2 array (k=1)'
101	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
102	ENDIF
103
104
105	! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
106	! ================================ whatever its orientation is)
107	! (USER: modify ahm3 and ahm4 following your desiderata)
108	! Here: ahm is proportional to the cube of the maximum of the gridspacing
109	! in the to horizontal direction
110
111	IF( ln_dynldf_bilap ) THEN
112
113	zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
114	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
115
116	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
117	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
118
119	za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
120	DO jj = 1, jpj
121	DO ji = 1, jpi
122	zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
123	zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
124	ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
125	ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
126	END DO
127	END DO
128
129	zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0
130	IF( ln_zco ) THEN ! z-coordinate, same profile everywhere
131	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
132	DO jk = 1, jpk
133	IF( fsdept(1,1,jk) <= zh ) THEN
134	zcoef(jk) = 1.e0
135	ELSE
136	zcoef(jk) = 1.e0 + ( zr - 1.e0 ) &
137	& * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) &
138	& / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) )
139	ENDIF
140	ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
141	ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
142	IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
143	END DO
144	ELSE ! partial steps or s-ccordinate
145	zc = MAXVAL( fsdept(:,:,jpkm1) )
146	IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain
147
148	zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) )
149	DO jk = 2, jpkm1
150	DO jj = 1, jpj
151	DO ji = 1, jpi
152	IF( fsdept(ji,jj,jk) <= zh ) THEN
153	ahm3(ji,jj,jk) = ahm3(ji,jj,1)
154	ahm4(ji,jj,jk) = ahm4(ji,jj,1)
155	ELSE
156	zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc
157	ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
158	ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
159	ENDIF
160	END DO
161	END DO
162	END DO
163	ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
164	ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
165	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
166	DO jk = 1, jpk
167	IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
168	END DO
169	ENDIF
170
171	! Control print
172	IF( lwp .AND. ld_print ) THEN
173	WRITE(numout,*)
174	WRITE(numout,*) 'inildf: ahm3 array at level 1'
175	CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
176	WRITE(numout,*)
177	WRITE(numout,*) 'inildf: ahm4 array at level 1'
178	CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
179	ENDIF
180	ENDIF
181	!
182	IF(lwp .AND. lflush) CALL flush(numout)
183	!
184	CALL wrk_dealloc( jpk, zcoef )
185	!
186	END SUBROUTINE ldf_dyn_c3d
187
188
189	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
190	!!----------------------------------------------------------------------
191	!! * ROUTINE ldf_dyn_c3d *
192	!!
193	!! ** Purpose : ORCA R1, R2 and R4 only
194	!!
195	!! ** Method : blah blah blah ....
196	!!----------------------------------------------------------------------
197	USE ldftra_oce, ONLY: aht0
198	USE iom
199	!!
200	LOGICAL, INTENT(in) :: ld_print ! If true, output arrays on numout
201	!!
202	INTEGER :: ji, jj, jk, jn ! dummy loop indices
203	INTEGER :: ii0, ii1, ij0, ij1 ! local integers
204	INTEGER :: inum, iim, ijm !
205	INTEGER :: ifreq, il1, il2, ij, ii
206	REAL(wp) :: zahmeq, zcoff, zcoft, zmsk ! local scalars
207	REAL(wp) :: zemax , zemin, zeref, zahmm
208	CHARACTER (len=15) :: clexp
209	INTEGER , POINTER, DIMENSION(:,:) :: icof
210	REAL(wp), POINTER, DIMENSION(: ) :: zcoef
211	REAL(wp), POINTER, DIMENSION(:,:) :: zahm0
212	!
213	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d ! temporary array to read ahmcoef file
214	LOGICAL :: tempmask( jpi,jpj) ! Temporary mask to avoid Cray compiler bug at cce 8.3.4
215	!!----------------------------------------------------------------------
216	!
217	CALL wrk_alloc( jpi , jpj , icof )
218	CALL wrk_alloc( jpk , zcoef )
219	CALL wrk_alloc( jpi , jpj , zahm0 )
220	!
221	IF(lwp) WRITE(numout,*)
222	IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
223	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
224	IF(lwp) WRITE(numout,*) ' orca R1, R2 or R4 configuration: reduced in the surface Eq. strip '
225
226	! Read 2d integer array to specify western boundary increase in the
227	! ===================== equatorial strip (20N-20S) defined at t-points
228	ALLOCATE( ztemp2d(jpi,jpj) )
229	ztemp2d(:,:) = 0.
230	CALL iom_open ( 'ahmcoef.nc', inum )
231	CALL iom_get ( inum, jpdom_data, 'icof', ztemp2d)
232	icof(:,:) = NINT(ztemp2d(:,:))
233	CALL iom_close( inum )
234	DEALLOCATE(ztemp2d)
235
236	! Set ahm1 and ahm2
237	! =================
238
239	! define ahm1 and ahm2 at the right grid point position
240	! (USER: modify ahm1 and ahm2 following your desiderata)
241	! biharmonic : ahm1 (ahm2) defined at u- (v-) point
242	! harmonic : ahm1 (ahm2) defined at t- (f-) point
243
244	! first level : as for 2D coefficients
245
246	! Decrease ahm to zahmeq m2/s in the tropics
247	! (from 90 to 20 degre: ahm = constant
248	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
249	! from 2.5 to 0 degre: ahm = constant
250	! symmetric in the south hemisphere)
251
252	IF( jp_cfg == 4 ) THEN
253	zahmeq = 5.0 * aht0
254	zahmm = min( 160000.0, ahm0)
255	zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
256	zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
257	tempmask(:,:) = .FALSE.
258	! Pre calculate mask for zeref since embedding the following
259	! term in the MAXVAL operation offends the Cray compiler for no
260	! justifiable reason under certain conditions.
261	tempmask(:,:) = (tmask(:,:,1) .GE. 0.5) .AND. (ABS(gphit(:,:)) .GT. 50.)
262	zeref = MAXVAL ( e1t(:,:) * e2t(:,:), tempmask(:,:) )
263
264	DO jj = 1, jpj
265	DO ji = 1, jpi
266	zmsk = e1t(ji,jj) * e2t(ji,jj)
267	IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
268	zahm0(ji,jj) = ahm0
269	ELSE
270	IF( zmsk .GE. zeref ) THEN
271	zahm0(ji,jj) = ahm0
272	ELSE
273	zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - &
274	& cos((rpi0.5(zmsk-zemin)/(zeref-zemin))))
275	ENDIF
276	ENDIF
277	END DO
278	END DO
279	ENDIF
280
281	IF( jp_cfg == 2 ) THEN
282	zahmeq = aht0
283	zahmm = ahm0
284	zahm0(:,:) = ahm0
285	ENDIF
286
287	IF( jp_cfg == 1 ) THEN
288	zahmeq = aht0 ! reduced to aht0 on equator; set to ahm0 if no tropical reduction is required
289	zahmm = ahm0
290	zahm0(:,:) = ahm0
291	ENDIF
292
293	DO jj = 1, jpj
294	DO ji = 1, jpi
295	IF( ABS(gphif(ji,jj)) >= 20.) THEN
296	ahm2(ji,jj,1) = zahm0(ji,jj)
297	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
298	ahm2(ji,jj,1) = zahmeq
299	ELSE
300	ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
301	& (1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
302	ENDIF
303	IF( ABS(gphit(ji,jj)) >= 20.) THEN
304	ahm1(ji,jj,1) = zahm0(ji,jj)
305	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
306	ahm1(ji,jj,1) = zahmeq
307	ELSE
308	ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
309	& (1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
310	ENDIF
311	END DO
312	END DO
313
314	! increase along western boundaries of equatorial strip
315	! t-point
316	DO jj = 1, jpjm1
317	DO ji = 1, jpim1
318	zcoft = float( icof(ji,jj) ) / 100.
319	ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
320	END DO
321	END DO
322	! f-point
323	icof(:,:) = icof(:,:) * tmask(:,:,1)
324	DO jj = 1, jpjm1
325	DO ji = 1, jpim1 ! NO vector opt.
326	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
327	IF( zmsk == 0. ) THEN
328	zcoff = 1.
329	ELSE
330	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
331	/ (zmsk * 100.)
332	ENDIF
333	ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
334	END DO
335	END DO
336
337	! other level: re-increase the coef in the deep ocean
338	!==================================================================
339	! Prior to v3.3, zcoeff was hardwired according to k-index jk.
340	!
341	! From v3.3 onwards this has been generalised to a function of
342	! depth so that it can be used with any number of levels.
343	!
344	! The function has been chosen to match the original values (shown
345	! in the following comments) when using the standard 31 ORCA levels.
346	! DO jk = 1, 21
347	! zcoef(jk) = 1._wp
348	! END DO
349	! zcoef(22) = 2._wp
350	! zcoef(23) = 3._wp
351	! zcoef(24) = 5._wp
352	! zcoef(25) = 7._wp
353	! zcoef(26) = 9._wp
354	! DO jk = 27, jpk
355	! zcoef(jk) = 10._wp
356	! END DO
357	!==================================================================
358
359	IF(lwp) THEN
360	WRITE(numout,*)
361	WRITE(numout,*) ' 1D zcoef array '
362	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
363	WRITE(numout,*)
364	WRITE(numout,*) ' jk zcoef '
365	ENDIF
366
367	DO jk=1, jpk
368	zcoef(jk) = 1.0_wp + NINT(9.0_wp*(gdept_1d(jk)-800.0_wp)/(3000.0_wp-800.0_wp))
369	zcoef(jk) = MIN(10.0_wp, MAX(1.0_wp, zcoef(jk)))
370	IF(lwp) WRITE(numout,'(4x,i3,6x,f7.3)') jk,zcoef(jk)
371	END DO
372
373	DO jk = 2, jpk
374	ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
375	ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
376	END DO
377
378	IF( jp_cfg == 4 ) THEN ! Limit AHM in Gibraltar strait
379	ij0 = 50 ; ij1 = 53
380	ii0 = 69 ; ii1 = 71
381	DO jk = 1, jpk
382	ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
383	ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
384	END DO
385	ENDIF
386	CALL lbc_lnk( ahm1, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)
387	CALL lbc_lnk( ahm2, 'F', 1. )
388
389
390	IF(lwp) THEN ! Control print
391	WRITE(numout,*)
392	WRITE(numout,*) ' 3D ahm1 array (k=1)'
393	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
394	WRITE(numout,*)
395	WRITE(numout,*) ' 3D ahm2 array (k=1)'
396	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
397	WRITE(numout,*)
398	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
399	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
400	ENDIF
401
402
403	! Set ahm3 and ahm4
404	! =================
405
406	! define ahm3 and ahm4 at the right grid point position
407	! initialization to a constant value
408	! (USER: modify ahm3 and ahm4 following your desiderata)
409	! harmonic isopycnal or geopotential:
410	! ahm3 (ahm4) defined at u- (v-) point
411	DO jk = 1, jpk
412	DO jj = 2, jpj
413	DO ji = 2, jpi
414	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
415	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
416	END DO
417	END DO
418	END DO
419	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
420	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
421
422	CALL lbc_lnk( ahm3, 'U', 1. ) ! Lateral boundary conditions (unchanged sign)
423	CALL lbc_lnk( ahm4, 'V', 1. )
424
425	! Control print
426
427	IF( lwp .AND. ld_print ) THEN
428	WRITE(numout,*)
429	WRITE(numout,*) ' ahm3 array level 1'
430	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
431	WRITE(numout,*)
432	WRITE(numout,*) ' ahm4 array level 1'
433	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
434	ENDIF
435	!
436	IF(lwp .AND. lflush) CALL flush(numout)
437	!
438	CALL wrk_dealloc( jpi , jpj , icof )
439	CALL wrk_dealloc( jpk , zcoef )
440	CALL wrk_dealloc( jpi , jpj , zahm0 )
441	!
442	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 11101

Download in other formats: