Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ldfdyn_c3d.h90 @ 6032

Last change on this file since 6032 was 6032, checked in by acosce, 2 years ago

create new configuration for IPSLCM5A2CHT model :

1- this configuration will allow to done paleo experiment (need to be done)
2- compilation will be done with script like v6.2 configurations

File size: 17.9 KB

Line
1	!!----------------------------------------------------------------------
2	!! * ldfdyn_c3d.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
7	!! $Id: ldfdyn_c3d.h90 5400 2015-06-10 15:29:08Z cbricaud $
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	! Add PALEORCA2 configuration -- JBL 08.02.2017
85	! ==============================================
86	IF( ( cp_cfg == "orca" .OR. cp_cfg == "paleorca" ) .AND. ( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
87	IF(lwp) WRITE(numout,*)
88	IF(lwp) WRITE(numout,*) ' ORCA R1, R2 or R4: overwrite the previous definition of ahm'
89	IF(lwp) WRITE(numout,*) ' ================='
90	CALL ldf_dyn_c3d_orca( ld_print )
91	ENDIF
92
93	ENDIF
94
95	! Control print
96	IF(lwp .AND. ld_print ) THEN
97	WRITE(numout,*)
98	WRITE(numout,*) ' 3D ahm1 array (k=1)'
99	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
100	WRITE(numout,*)
101	WRITE(numout,*) ' 3D ahm2 array (k=1)'
102	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
103	ENDIF
104
105
106	! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
107	! ================================ whatever its orientation is)
108	! (USER: modify ahm3 and ahm4 following your desiderata)
109	! Here: ahm is proportional to the cube of the maximum of the gridspacing
110	! in the to horizontal direction
111
112	IF( ln_dynldf_bilap ) THEN
113
114	zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
115	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
116
117	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
118	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
119
120	za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
121	DO jj = 1, jpj
122	DO ji = 1, jpi
123	zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
124	zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
125	ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
126	ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
127	END DO
128	END DO
129
130	zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0
131	IF( ln_zco ) THEN ! z-coordinate, same profile everywhere
132	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
133	DO jk = 1, jpk
134	IF( fsdept(1,1,jk) <= zh ) THEN
135	zcoef(jk) = 1.e0
136	ELSE
137	zcoef(jk) = 1.e0 + ( zr - 1.e0 ) &
138	& * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) &
139	& / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) )
140	ENDIF
141	ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
142	ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
143	IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
144	END DO
145	ELSE ! partial steps or s-ccordinate
146	zc = MAXVAL( fsdept(:,:,jpkm1) )
147	IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain
148
149	zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) )
150	DO jk = 2, jpkm1
151	DO jj = 1, jpj
152	DO ji = 1, jpi
153	IF( fsdept(ji,jj,jk) <= zh ) THEN
154	ahm3(ji,jj,jk) = ahm3(ji,jj,1)
155	ahm4(ji,jj,jk) = ahm4(ji,jj,1)
156	ELSE
157	zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc
158	ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
159	ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
160	ENDIF
161	END DO
162	END DO
163	END DO
164	ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
165	ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
166	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
167	DO jk = 1, jpk
168	IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
169	END DO
170	ENDIF
171
172	! Control print
173	IF( lwp .AND. ld_print ) THEN
174	WRITE(numout,*)
175	WRITE(numout,*) 'inildf: ahm3 array at level 1'
176	CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
177	WRITE(numout,*)
178	WRITE(numout,*) 'inildf: ahm4 array at level 1'
179	CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
180	ENDIF
181	ENDIF
182	!
183	CALL wrk_dealloc( jpk, zcoef )
184	!
185	END SUBROUTINE ldf_dyn_c3d
186
187
188	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
189	!!----------------------------------------------------------------------
190	!! * ROUTINE ldf_dyn_c3d *
191	!!
192	!! ** Purpose : ORCA R1, R2 and R4 only
193	!!
194	!! ** Method : blah blah blah ....
195	!!----------------------------------------------------------------------
196	USE ldftra_oce, ONLY: aht0
197	USE iom
198	!!
199	LOGICAL, INTENT(in) :: ld_print ! If true, output arrays on numout
200	!!
201	INTEGER :: ji, jj, jk, jn ! dummy loop indices
202	INTEGER :: ii0, ii1, ij0, ij1 ! local integers
203	INTEGER :: inum, iim, ijm !
204	INTEGER :: ifreq, il1, il2, ij, ii
205	REAL(wp) :: zahmeq, zcoff, zcoft, zmsk ! local scalars
206	REAL(wp) :: zemax , zemin, zeref, zahmm
207	CHARACTER (len=15) :: clexp
208	INTEGER , POINTER, DIMENSION(:,:) :: icof
209	REAL(wp), POINTER, DIMENSION(: ) :: zcoef
210	REAL(wp), POINTER, DIMENSION(:,:) :: zahm0
211	!
212	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d ! temporary array to read ahmcoef file
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	zeref = MAXVAL ( e1t(:,:) * e2t(:,:), &
256	& tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
257
258	DO jj = 1, jpj
259	DO ji = 1, jpi
260	zmsk = e1t(ji,jj) * e2t(ji,jj)
261	IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
262	zahm0(ji,jj) = ahm0
263	ELSE
264	IF( zmsk .GE. zeref ) THEN
265	zahm0(ji,jj) = ahm0
266	ELSE
267	zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - &
268	& cos((rpi0.5(zmsk-zemin)/(zeref-zemin))))
269	ENDIF
270	ENDIF
271	END DO
272	END DO
273	ENDIF
274
275	IF( jp_cfg == 2 ) THEN
276	zahmeq = aht0
277	zahmm = ahm0
278	zahm0(:,:) = ahm0
279	ENDIF
280
281	IF( jp_cfg == 1 ) THEN
282	zahmeq = aht0 ! reduced to aht0 on equator; set to ahm0 if no tropical reduction is required
283	zahmm = ahm0
284	zahm0(:,:) = ahm0
285	ENDIF
286
287	DO jj = 1, jpj
288	DO ji = 1, jpi
289	IF( ABS(gphif(ji,jj)) >= 20.) THEN
290	ahm2(ji,jj,1) = zahm0(ji,jj)
291	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
292	ahm2(ji,jj,1) = zahmeq
293	ELSE
294	ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
295	& (1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
296	ENDIF
297	IF( ABS(gphit(ji,jj)) >= 20.) THEN
298	ahm1(ji,jj,1) = zahm0(ji,jj)
299	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
300	ahm1(ji,jj,1) = zahmeq
301	ELSE
302	ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
303	& (1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
304	ENDIF
305	END DO
306	END DO
307
308	! increase along western boundaries of equatorial strip
309	! t-point
310	DO jj = 1, jpjm1
311	DO ji = 1, jpim1
312	zcoft = float( icof(ji,jj) ) / 100.
313	ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
314	END DO
315	END DO
316	! f-point
317	icof(:,:) = icof(:,:) * tmask(:,:,1)
318	DO jj = 1, jpjm1
319	DO ji = 1, jpim1 ! NO vector opt.
320	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
321	IF( zmsk == 0. ) THEN
322	zcoff = 1.
323	ELSE
324	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
325	/ (zmsk * 100.)
326	ENDIF
327	ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
328	END DO
329	END DO
330
331	! other level: re-increase the coef in the deep ocean
332	!==================================================================
333	! Prior to v3.3, zcoeff was hardwired according to k-index jk.
334	!
335	! From v3.3 onwards this has been generalised to a function of
336	! depth so that it can be used with any number of levels.
337	!
338	! The function has been chosen to match the original values (shown
339	! in the following comments) when using the standard 31 ORCA levels.
340	! DO jk = 1, 21
341	! zcoef(jk) = 1._wp
342	! END DO
343	! zcoef(22) = 2._wp
344	! zcoef(23) = 3._wp
345	! zcoef(24) = 5._wp
346	! zcoef(25) = 7._wp
347	! zcoef(26) = 9._wp
348	! DO jk = 27, jpk
349	! zcoef(jk) = 10._wp
350	! END DO
351	!==================================================================
352
353	IF(lwp) THEN
354	WRITE(numout,*)
355	WRITE(numout,*) ' 1D zcoef array '
356	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
357	WRITE(numout,*)
358	WRITE(numout,*) ' jk zcoef '
359	ENDIF
360
361	DO jk=1, jpk
362	zcoef(jk) = 1.0_wp + NINT(9.0_wp*(gdept_1d(jk)-800.0_wp)/(3000.0_wp-800.0_wp))
363	zcoef(jk) = MIN(10.0_wp, MAX(1.0_wp, zcoef(jk)))
364	IF(lwp) WRITE(numout,'(4x,i3,6x,f7.3)') jk,zcoef(jk)
365	END DO
366
367	DO jk = 2, jpk
368	ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
369	ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
370	END DO
371
372	IF( jp_cfg == 4 ) THEN ! Limit AHM in Gibraltar strait
373	ij0 = 50 ; ij1 = 53
374	ii0 = 69 ; ii1 = 71
375	DO jk = 1, jpk
376	ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
377	ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
378	END DO
379	ENDIF
380	CALL lbc_lnk( ahm1, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)
381	CALL lbc_lnk( ahm2, 'F', 1. )
382
383
384	IF(lwp) THEN ! Control print
385	WRITE(numout,*)
386	WRITE(numout,*) ' 3D ahm1 array (k=1)'
387	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
388	WRITE(numout,*)
389	WRITE(numout,*) ' 3D ahm2 array (k=1)'
390	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
391	WRITE(numout,*)
392	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
393	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
394	ENDIF
395
396
397	! Set ahm3 and ahm4
398	! =================
399
400	! define ahm3 and ahm4 at the right grid point position
401	! initialization to a constant value
402	! (USER: modify ahm3 and ahm4 following your desiderata)
403	! harmonic isopycnal or geopotential:
404	! ahm3 (ahm4) defined at u- (v-) point
405	DO jk = 1, jpk
406	DO jj = 2, jpj
407	DO ji = 2, jpi
408	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
409	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
410	END DO
411	END DO
412	END DO
413	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
414	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
415
416	CALL lbc_lnk( ahm3, 'U', 1. ) ! Lateral boundary conditions (unchanged sign)
417	CALL lbc_lnk( ahm4, 'V', 1. )
418
419	! Control print
420
421	IF( lwp .AND. ld_print ) THEN
422	WRITE(numout,*)
423	WRITE(numout,*) ' ahm3 array level 1'
424	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
425	WRITE(numout,*)
426	WRITE(numout,*) ' ahm4 array level 1'
427	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
428	ENDIF
429	!
430	CALL wrk_dealloc( jpi , jpj , icof )
431	CALL wrk_dealloc( jpk , zcoef )
432	CALL wrk_dealloc( jpi , jpj , zahm0 )
433	!
434	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format