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_c3d.h90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

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

Last change on this file since 9975 was 9975, checked in by marc, 6 years ago
GMED ticket 400. Fixes to make the code more portable.
File size: 18.2 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	CALL wrk_dealloc( jpk, zcoef )
183	!
184	END SUBROUTINE ldf_dyn_c3d
185
186
187	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
188	!!----------------------------------------------------------------------
189	!! * ROUTINE ldf_dyn_c3d *
190	!!
191	!! ** Purpose : ORCA R1, R2 and R4 only
192	!!
193	!! ** Method : blah blah blah ....
194	!!----------------------------------------------------------------------
195	USE ldftra_oce, ONLY: aht0
196	USE iom
197	!!
198	LOGICAL, INTENT(in) :: ld_print ! If true, output arrays on numout
199	!!
200	INTEGER :: ji, jj, jk, jn ! dummy loop indices
201	INTEGER :: ii0, ii1, ij0, ij1 ! local integers
202	INTEGER :: inum, iim, ijm !
203	INTEGER :: ifreq, il1, il2, ij, ii
204	REAL(wp) :: zahmeq, zcoff, zcoft, zmsk ! local scalars
205	REAL(wp) :: zemax , zemin, zeref, zahmm
206	CHARACTER (len=15) :: clexp
207	INTEGER , POINTER, DIMENSION(:,:) :: icof
208	REAL(wp), POINTER, DIMENSION(: ) :: zcoef
209	REAL(wp), POINTER, DIMENSION(:,:) :: zahm0
210	!
211	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d ! temporary array to read ahmcoef file
212	LOGICAL :: tempmask( jpi,jpj) ! Temporary mask to avoid Cray compiler bug at cce 8.3.4
213	!!----------------------------------------------------------------------
214	!
215	CALL wrk_alloc( jpi , jpj , icof )
216	CALL wrk_alloc( jpk , zcoef )
217	CALL wrk_alloc( jpi , jpj , zahm0 )
218	!
219	IF(lwp) WRITE(numout,*)
220	IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
221	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
222	IF(lwp) WRITE(numout,*) ' orca R1, R2 or R4 configuration: reduced in the surface Eq. strip '
223
224	! Read 2d integer array to specify western boundary increase in the
225	! ===================== equatorial strip (20N-20S) defined at t-points
226	ALLOCATE( ztemp2d(jpi,jpj) )
227	ztemp2d(:,:) = 0.
228	CALL iom_open ( 'ahmcoef.nc', inum )
229	CALL iom_get ( inum, jpdom_data, 'icof', ztemp2d)
230	icof(:,:) = NINT(ztemp2d(:,:))
231	CALL iom_close( inum )
232	DEALLOCATE(ztemp2d)
233
234	! Set ahm1 and ahm2
235	! =================
236
237	! define ahm1 and ahm2 at the right grid point position
238	! (USER: modify ahm1 and ahm2 following your desiderata)
239	! biharmonic : ahm1 (ahm2) defined at u- (v-) point
240	! harmonic : ahm1 (ahm2) defined at t- (f-) point
241
242	! first level : as for 2D coefficients
243
244	! Decrease ahm to zahmeq m2/s in the tropics
245	! (from 90 to 20 degre: ahm = constant
246	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
247	! from 2.5 to 0 degre: ahm = constant
248	! symmetric in the south hemisphere)
249
250	IF( jp_cfg == 4 ) THEN
251	zahmeq = 5.0 * aht0
252	zahmm = min( 160000.0, ahm0)
253	zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
254	zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
255	tempmask(:,:) = .FALSE.
256	! Pre calculate mask for zeref since embedding the following
257	! term in the MAXVAL operation offends the Cray compiler for no
258	! justifiable reason under certain conditions.
259	tempmask(:,:) = (tmask(:,:,1) .GE. 0.5) .AND. (ABS(gphit(:,:)) .GT. 50.)
260	zeref = MAXVAL ( e1t(:,:) * e2t(:,:), tempmask(:,:) )
261
262	DO jj = 1, jpj
263	DO ji = 1, jpi
264	zmsk = e1t(ji,jj) * e2t(ji,jj)
265	IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
266	zahm0(ji,jj) = ahm0
267	ELSE
268	IF( zmsk .GE. zeref ) THEN
269	zahm0(ji,jj) = ahm0
270	ELSE
271	zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - &
272	& cos((rpi0.5(zmsk-zemin)/(zeref-zemin))))
273	ENDIF
274	ENDIF
275	END DO
276	END DO
277	ENDIF
278
279	IF( jp_cfg == 2 ) THEN
280	zahmeq = aht0
281	zahmm = ahm0
282	zahm0(:,:) = ahm0
283	ENDIF
284
285	IF( jp_cfg == 1 ) THEN
286	zahmeq = aht0 ! reduced to aht0 on equator; set to ahm0 if no tropical reduction is required
287	zahmm = ahm0
288	zahm0(:,:) = ahm0
289	ENDIF
290
291	DO jj = 1, jpj
292	DO ji = 1, jpi
293	IF( ABS(gphif(ji,jj)) >= 20.) THEN
294	ahm2(ji,jj,1) = zahm0(ji,jj)
295	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
296	ahm2(ji,jj,1) = zahmeq
297	ELSE
298	ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
299	& (1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
300	ENDIF
301	IF( ABS(gphit(ji,jj)) >= 20.) THEN
302	ahm1(ji,jj,1) = zahm0(ji,jj)
303	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
304	ahm1(ji,jj,1) = zahmeq
305	ELSE
306	ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
307	& (1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
308	ENDIF
309	END DO
310	END DO
311
312	! increase along western boundaries of equatorial strip
313	! t-point
314	DO jj = 1, jpjm1
315	DO ji = 1, jpim1
316	zcoft = float( icof(ji,jj) ) / 100.
317	ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
318	END DO
319	END DO
320	! f-point
321	icof(:,:) = icof(:,:) * tmask(:,:,1)
322	DO jj = 1, jpjm1
323	DO ji = 1, jpim1 ! NO vector opt.
324	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
325	IF( zmsk == 0. ) THEN
326	zcoff = 1.
327	ELSE
328	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
329	/ (zmsk * 100.)
330	ENDIF
331	ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
332	END DO
333	END DO
334
335	! other level: re-increase the coef in the deep ocean
336	!==================================================================
337	! Prior to v3.3, zcoeff was hardwired according to k-index jk.
338	!
339	! From v3.3 onwards this has been generalised to a function of
340	! depth so that it can be used with any number of levels.
341	!
342	! The function has been chosen to match the original values (shown
343	! in the following comments) when using the standard 31 ORCA levels.
344	! DO jk = 1, 21
345	! zcoef(jk) = 1._wp
346	! END DO
347	! zcoef(22) = 2._wp
348	! zcoef(23) = 3._wp
349	! zcoef(24) = 5._wp
350	! zcoef(25) = 7._wp
351	! zcoef(26) = 9._wp
352	! DO jk = 27, jpk
353	! zcoef(jk) = 10._wp
354	! END DO
355	!==================================================================
356
357	IF(lwp) THEN
358	WRITE(numout,*)
359	WRITE(numout,*) ' 1D zcoef array '
360	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
361	WRITE(numout,*)
362	WRITE(numout,*) ' jk zcoef '
363	ENDIF
364
365	DO jk=1, jpk
366	zcoef(jk) = 1.0_wp + NINT(9.0_wp*(gdept_1d(jk)-800.0_wp)/(3000.0_wp-800.0_wp))
367	zcoef(jk) = MIN(10.0_wp, MAX(1.0_wp, zcoef(jk)))
368	IF(lwp) WRITE(numout,'(4x,i3,6x,f7.3)') jk,zcoef(jk)
369	END DO
370
371	DO jk = 2, jpk
372	ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
373	ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
374	END DO
375
376	IF( jp_cfg == 4 ) THEN ! Limit AHM in Gibraltar strait
377	ij0 = 50 ; ij1 = 53
378	ii0 = 69 ; ii1 = 71
379	DO jk = 1, jpk
380	ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
381	ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
382	END DO
383	ENDIF
384	CALL lbc_lnk( ahm1, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)
385	CALL lbc_lnk( ahm2, 'F', 1. )
386
387
388	IF(lwp) THEN ! Control print
389	WRITE(numout,*)
390	WRITE(numout,*) ' 3D ahm1 array (k=1)'
391	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
392	WRITE(numout,*)
393	WRITE(numout,*) ' 3D ahm2 array (k=1)'
394	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
395	WRITE(numout,*)
396	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
397	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
398	ENDIF
399
400
401	! Set ahm3 and ahm4
402	! =================
403
404	! define ahm3 and ahm4 at the right grid point position
405	! initialization to a constant value
406	! (USER: modify ahm3 and ahm4 following your desiderata)
407	! harmonic isopycnal or geopotential:
408	! ahm3 (ahm4) defined at u- (v-) point
409	DO jk = 1, jpk
410	DO jj = 2, jpj
411	DO ji = 2, jpi
412	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
413	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
414	END DO
415	END DO
416	END DO
417	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
418	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
419
420	CALL lbc_lnk( ahm3, 'U', 1. ) ! Lateral boundary conditions (unchanged sign)
421	CALL lbc_lnk( ahm4, 'V', 1. )
422
423	! Control print
424
425	IF( lwp .AND. ld_print ) THEN
426	WRITE(numout,*)
427	WRITE(numout,*) ' ahm3 array level 1'
428	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
429	WRITE(numout,*)
430	WRITE(numout,*) ' ahm4 array level 1'
431	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
432	ENDIF
433	!
434	CALL wrk_dealloc( jpi , jpj , icof )
435	CALL wrk_dealloc( jpk , zcoef )
436	CALL wrk_dealloc( jpi , jpj , zahm0 )
437	!
438	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format