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.

trcdmp.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: trunk/NEMO/TOP_SRC/TRP/trcdmp.F90 @ 281

Last change on this file since 281 was 274, checked in by opalod, 19 years ago
nemo_v1_update_005:RB: update headers for the TOP component.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 35.2 KB

Line
1	MODULE trcdmp
2	!!======================================================================
3	!! * MODULE trcdmp *
4	!! Ocean physics: internal restoring trend on passive tracers
5	!!======================================================================
6	#if defined key_passivetrc && defined key_trcdmp
7	!!----------------------------------------------------------------------
8	!! key_trcdmp internal damping
9	!!----------------------------------------------------------------------
10	!! trc_dmp : update the tracer trend with the internal damping
11	!! trc_dmp_init : initialization, namlist read, parameters control
12	!! trccof_zoom : restoring coefficient for zoom domain
13	!! trccof : restoring coefficient for global domain
14	!! cofdis : compute the distance to the coastline
15	!!----------------------------------------------------------------------
16	!! * Modules used
17	USE oce_trc ! ocean dynamics and tracers variables
18	USE trc ! ocean passive tracers variables
19	USE trctrp_lec ! passive tracers transport
20	USE trcdta
21
22	IMPLICIT NONE
23	PRIVATE
24
25	!! * Routine accessibility
26	PUBLIC trc_dmp ! routine called by step.F90
27
28	!! * Shared module variables
29	LOGICAL , PUBLIC, PARAMETER :: lk_trcdmp = .TRUE. !: internal damping flag
30
31	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: &
32	restotr ! restoring coeff. on tracers (s-1)
33
34	!! * Substitutions
35	# include "passivetrc_substitute.h90"
36	!!----------------------------------------------------------------------
37	!! TOP 1.0, LOCEAN-IPSL (2005)
38	!! $Header$
39	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
40	!!----------------------------------------------------------------------
41
42	CONTAINS
43
44	SUBROUTINE trc_dmp( kt )
45	!!----------------------------------------------------------------------
46	!! * ROUTINE trc_dmp *
47	!!
48	!! ** Purpose : Compute the passive tracer trend due to a newtonian damping
49	!! of the tracer field towards given data field and add it to the
50	!! general tracer trends.
51	!!
52	!! ** Method : Newtonian damping towards trdta computed
53	!! and add to the general tracer trends:
54	!! trn = tra + restotr * (trdta - trb)
55	!! The trend is computed either throughout the water column
56	!! (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
57	!! below the well mixed layer (nlmdmptr=2)
58	!!
59	!! ** Action : - update the tracer trends tra with the newtonian
60	!! damping trends.
61	!! - save the trends in trtrd ('key_trc_diatrd')
62	!!
63	!! History :
64	!! 7.0 ! (G. Madec) Original code
65	!! ! 96-01 (G. Madec)
66	!! ! 97-05 (H. Loukos) adapted for passive tracers
67	!! 8.5 ! 02-08 (G. Madec ) free form + modules
68	!! 9.0 ! 04-03 (C. Ethe) free form + modules
69	!!----------------------------------------------------------------------
70	!! * Arguments
71	INTEGER, INTENT( in ) :: kt ! ocean time-step index
72
73	!! * Local declarations
74	INTEGER :: ji, jj, jk, jn ! dummy loop indices
75	REAL(wp) :: ztest, ztra, zdt ! temporary scalars
76	!!----------------------------------------------------------------------
77
78	! 0. Initialization (first time-step only)
79	! --------------
80	IF( kt == nittrc000 ) CALL trc_dmp_init
81
82	! 1. Newtonian damping trends on tracer fields
83	! --------------------------------------------
84	! compute the newtonian damping trends depending on nmldmptr
85
86	!!! zdt = rdt * FLOAT( ndttrc )
87
88	! Initialize the input fields for newtonian damping
89	CALL dta_trc( kt )
90
91	DO jn = 1, jptra
92
93	IF( lutini(jn) ) THEN
94
95	SELECT CASE ( nmldmptr )
96
97	CASE( 0 ) ! newtonian damping throughout the water column
98
99	DO jk = 1, jpkm1
100	DO jj = 2, jpjm1
101	DO ji = fs_2, fs_jpim1 ! vector opt.
102	ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
103	! add the trends to the general tracer trends
104	!! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
105	tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
106	# if defined key_trc_diatrd
107	! save the trends for diagnostics
108	trtrd(ji,jj,jk,jn,jpdiatrc) = ztra
109	# endif
110	END DO
111	END DO
112	END DO
113
114	CASE ( 1 ) ! no damping in the turbocline (avt > 5 cm2/s)
115	DO jk = 1, jpkm1
116	DO jj = 2, jpjm1
117	DO ji = fs_2, fs_jpim1 ! vector opt.
118	ztest = avt(ji,jj,jk) - 5.e-4
119	IF( ztest < 0. ) THEN
120	ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
121	ELSE
122	ztra = 0.e0
123	ENDIF
124	! add the trends to the general tracer trends
125	!! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
126	tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
127	# if defined key_trc_diatrd
128	! save the trends for diagnostics
129	trtrd(ji,jj,jk,jn,jpdiatrc) = ztra
130	# endif
131	END DO
132	END DO
133	END DO
134
135	CASE ( 2 ) ! no damping in the mixed layer
136	DO jk = 1, jpkm1
137	DO jj = 2, jpjm1
138	DO ji = fs_2, fs_jpim1 ! vector opt.
139	IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
140	ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
141	ELSE
142	ztra = 0.e0
143	ENDIF
144	! add the trends to the general tracer trends
145	!! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
146	tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
147	# if defined key_trc_diatrd
148	! save the trends for diagnostics
149	trtrd(ji,jj,jk,jn,jpdiatrc) = ztra
150	# endif
151	END DO
152	END DO
153	END DO
154
155	END SELECT
156
157	ENDIF
158
159	IF(ln_ctl) THEN ! print mean trends (used for debugging)
160	ztra = SUM( tra(2:nictle,2:njctle,1:jpkm1,jn) * tmask(2:nictle,2:njctle,1:jpkm1) )
161	WRITE(numout,*) ' trc/dmp - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn)
162	tra_ctl(jn) = ztra
163	ENDIF
164
165	END DO
166
167	trb(:,:,:,:) = trn(:,:,:,:)
168
169	END SUBROUTINE trc_dmp
170
171
172	SUBROUTINE trc_dmp_init
173	!!----------------------------------------------------------------------
174	!! * ROUTINE trc_dmp_init *
175	!!
176	!! ** Purpose : Initialization for the newtonian damping
177	!!
178	!! ** Method : read the nammbf namelist and check the parameters
179	!! called by trc_dmp at the first timestep (nit000)
180	!!
181	!! History :
182	!! 8.5 ! 02-08 (G. Madec) Original code
183	!!----------------------------------------------------------------------
184
185	SELECT CASE ( ndmptr )
186
187	CASE ( -1 ) ! ORCA: damping in Red & Med Seas only
188	IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only'
189
190	CASE ( 1:90 ) ! Damping poleward of 'ndmptr' degrees
191	IF(lwp) WRITE(numout,*) ' tracer damping poleward of', ndmptr, ' degrees'
192
193	CASE DEFAULT
194	IF(lwp) WRITE(numout,cform_err)
195	IF(lwp) WRITE(numout,*) ' bad flag value for ndmptr = ', ndmptr
196	nstop = nstop + 1
197
198	END SELECT
199
200
201	SELECT CASE ( nmldmptr )
202
203	CASE ( 0 ) ! newtonian damping throughout the water column
204	IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column'
205
206	CASE ( 1 ) ! no damping in the turbocline (avt > 5 cm2/s)
207	IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline'
208
209	CASE ( 2 ) ! no damping in the mixed layer
210	IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer'
211
212	CASE DEFAULT
213	IF(lwp) WRITE(numout,cform_err)
214	IF(lwp) WRITE(numout,*) ' bad flag value for nmldmptr = ', nmldmptr
215	nstop = nstop + 1
216
217	END SELECT
218
219
220	! 3. Damping coefficients initialization
221	! --------------------------------------
222
223	IF( lzoom ) THEN
224	CALL trccof_zoom
225	ELSE
226	CALL trccof
227	ENDIF
228
229	END SUBROUTINE trc_dmp_init
230
231
232	SUBROUTINE trccof_zoom
233	!!----------------------------------------------------------------------
234	!! * ROUTINE trccof_zoom *
235	!!
236	!! ** Purpose : Compute the damping coefficient for zoom domain
237	!!
238	!! ** Method : - set along closed boundary due to zoom a damping over
239	!! 6 points with a max time scale of 5 days.
240	!! - ORCA arctic/antarctic zoom: set the damping along
241	!! south/north boundary over a latitude strip.
242	!!
243	!! ** Action : - restotr, the damping coeff. passive tracers
244	!!
245	!! History :
246	!! 9.0 ! 03-09 (G. Madec) Original code
247	!! 9.0 ! 04-03 (C. Ethe) adapted for passive tracers
248	!!----------------------------------------------------------------------
249	!! * Local declarations
250	INTEGER :: ji, jj, jk, jn ! dummy loop indices
251	REAL(wp) :: &
252	zlat, zlat0, zlat1, zlat2 ! temporary scalar
253	REAL(wp), DIMENSION(6) :: &
254	zfact ! temporary workspace
255	!!----------------------------------------------------------------------
256
257	zfact(1) = 1.
258	zfact(2) = 1.
259	zfact(3) = 11./12.
260	zfact(4) = 8./12.
261	zfact(5) = 4./12.
262	zfact(6) = 1./12.
263	zfact(:) = zfact(:) / ( 5. * rday ) ! 5 days max restoring time scale
264
265	restotr(:,:,:,:) = 0.e0
266
267	! damping along the forced closed boundary over 6 grid-points
268	DO jn = 1, 6
269	IF( lzoom_w ) restotr( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : , : ) = zfact(jn) ! west closed
270	IF( lzoom_s ) restotr( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : , : ) = zfact(jn) ! south closed
271	IF( lzoom_e ) restotr( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : , : ) &
272	& = zfact(jn) ! east closed
273	IF( lzoom_n ) restotr( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : , : ) &
274	& = zfact(jn) ! north closed
275	END DO
276
277
278	IF( lzoom_arct .AND. lzoom_anta ) THEN
279
280	! ====================================================
281	! ORCA configuration : arctic zoom or antarctic zoom
282	! ====================================================
283
284	IF(lwp) WRITE(numout,*)
285	IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' trccof_zoom : ORCA Arctic zoom'
286	IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' trccof_zoom : ORCA Antarctic zoom'
287	IF(lwp) WRITE(numout,*)
288
289	! ... Initialization :
290	! zlat0 : latitude strip where resto decreases
291	! zlat1 : resto = 1 before zlat1
292	! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
293	restotr(:,:,:,:) = 0.e0
294	zlat0 = 10.
295	zlat1 = 30.
296	zlat2 = zlat1 + zlat0
297
298	! ... Compute arrays resto ; value for internal damping : 5 days
299	DO jn = 1, jptra
300	DO jk = 2, jpkm1
301	DO jj = 1, jpj
302	DO ji = 1, jpi
303	zlat = ABS( gphit(ji,jj) )
304	IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
305	restotr(ji,jj,jk,jn) = 0.5 * ( 1./(5.rday) ) &
306	( 1. - COS(rpi*(zlat2-zlat)/zlat0) )
307	ELSE IF ( zlat < zlat1 ) THEN
308	restotr(ji,jj,jk,jn) = 1./(5.*rday)
309	ENDIF
310	END DO
311	END DO
312	END DO
313	END DO
314
315	ENDIF
316
317	! ... Mask resto array
318	DO jn = 1, jptra
319	restotr(:,:,:,jn) = restotr(:,:,:,jn) * tmask(:,:,:)
320	END DO
321
322
323	END SUBROUTINE trccof_zoom
324
325	SUBROUTINE trccof
326	!!----------------------------------------------------------------------
327	!! * ROUTINE trccof *
328	!!
329	!! ** Purpose : Compute the damping coefficient
330	!!
331	!! ** Method : Arrays defining the damping are computed for each grid
332	!! point passive tracers (restotr)
333	!! Damping depends on distance to coast, depth and latitude
334	!!
335	!! ** Action : - restotr, the damping coeff. for passive tracers
336	!!
337	!! History :
338	!! 5.0 ! 91-03 (O. Marti, G. Madec) Original code
339	!! ! 92-06 (M. Imbard) doctor norme
340	!! ! 96-01 (G. Madec) statement function for e3
341	!! ! 98-07 (M. Imbard, G. Madec) ORCA version
342	!! ! 00-08 (G. Madec, D. Ludicone)
343	!! 8.2 ! 04-03 (H. Loukos) adapted for passive tracers
344	!! ! 04-02 (O. Aumont, C. Ethe) rewritten for debuging and update
345	!!----------------------------------------------------------------------
346	!! * Modules used
347	USE ioipsl
348
349	!! * Local declarations
350	INTEGER :: ji, jj, jk, je, jn ! dummy loop indices
351	INTEGER, PARAMETER :: jpmois=1
352	INTEGER :: ipi, ipj, ipk ! temporary integers
353	INTEGER :: ii0, ii1, ij0, ij1 ! " "
354	INTEGER :: &
355	idmp, & ! logical unit for file restoring damping term
356	icot ! logical unit for file distance to the coast
357	INTEGER :: itime, istep(jpmois), ie
358	LOGICAL :: llbon
359	CHARACTER (len=32) :: clname, clname2, clname3
360	REAL(wp) :: &
361	zdate0, zinfl, zlon, & ! temporary scalars
362	zlat, zlat0, zlat1, zlat2, & ! " "
363	zsdmp, zbdmp ! " "
364	REAL(wp), DIMENSION(jpk) :: &
365	zdept, zhfac
366	REAL(wp), DIMENSION(jpi,jpj) :: &
367	zmrs, zlamt, zphit
368	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
369	zdct
370	!!----------------------------------------------------------------------
371
372	! ====================================
373	! ORCA configuration : global domain
374	! ====================================
375
376	IF(lwp) WRITE(numout,*)
377	IF(lwp) WRITE(numout,*) ' trccof : Global domain of ORCA'
378	IF(lwp) WRITE(numout,*) ' ------------------------------'
379
380
381	! ... Initialization :
382	! zdct() : distant to the coastline
383	! resto() : array of restoring coeff.
384
385	zdct (:,:,:) = 0.e0
386	restotr(:,:,:,:) = 0.e0
387
388
389	IF ( ndmptr > 0 ) THEN
390
391	! ------------------------------------
392	! Damping poleward of 'ndmptr' degrees
393	! ------------------------------------
394
395	IF(lwp) WRITE(numout,*)
396	IF(lwp) WRITE(numout,*) ' Damping poleward of ', ndmptr,' deg.'
397	IF(lwp) WRITE(numout,*)
398
399	! ... Distance to coast (zdct)
400
401	! ... Test the existance of distance-to-coast file
402	itime = jpmois
403	ipi = jpiglo
404	ipj = jpjglo
405	ipk = jpk
406	clname = 'dist.coast.trc'
407	DO je = 1,32
408	IF( clname(je:je) == ' ' ) go to 140
409	END DO
410	140 CONTINUE
411	ie = je
412	clname2 = clname(1:ie-1)//".nc"
413	INQUIRE( FILE = clname2, EXIST = llbon )
414
415	IF ( llbon ) THEN
416
417	! ... Read file distance to coast if possible
418	CALL flinopen( clname, mig(1), nlci, mjg(1), nlcj, .FALSE., &
419	ipi, ipj, ipk, zlamt, zphit, zdept, jpmois, &
420	istep, zdate0, rdt, icot )
421	CALL flinget( icot, 'Tcoast', jpidta, jpjdta, jpk, &
422	jpmois, 1, 1, mig(1), nlci, mjg(1), nlcj, zdct(1:nlci,1:nlcj,1:jpk) )
423	CALL flinclo( icot )
424	IF(lwp)WRITE(numout,) ' * : File trc.dist.coast.nc read'
425
426	ELSE
427
428	! ... Compute and save the distance-to-coast array (output in zdct)
429	CALL cofdis ( zdct )
430
431	ENDIF
432
433	! ... Compute arrays resto
434	! zinfl : distance of influence for damping term
435	! zlat0 : latitude strip where resto decreases
436	! zlat1 : resto = 0 between -zlat1 and zlat1
437	! zlat2 : resto increases from 0 to 1 between \|zlat1\| and \|zlat2\|
438	! and resto = 1 between \|zlat2\| and 90 deg.
439	zinfl = 1000.e3
440	zlat0 = 10
441	zlat1 = ndmptr
442	zlat2 = zlat1 + zlat0
443
444	DO jn = 1, jptra
445	DO jj = 1, jpj
446	DO ji = 1, jpi
447	zlat = ABS( gphit(ji,jj) )
448	IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
449	restotr(ji,jj,1,jn) = 0.5 * ( 1. - COS(rpi*(zlat-zlat1)/zlat0 ) )
450	ELSEIF ( zlat > zlat2 ) THEN
451	restotr(ji,jj,1,jn) = 1.
452	ENDIF
453	END DO
454	END DO
455	END DO
456
457	! ... North Indian ocean (20N/30N x 45E/100E) : resto=0
458	IF ( ndmptr == 20 ) THEN
459	DO jn = 1, jptra
460	DO jj = 1, jpj
461	DO ji = 1, jpi
462	zlat = gphit(ji,jj)
463	zlon = MOD( glamt(ji,jj), 360. )
464	IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. &
465	45. < zlon .AND. zlon < 100. ) THEN
466	restotr(ji,jj,1,jn) = 0.
467	ENDIF
468	END DO
469	END DO
470	END DO
471	ENDIF
472
473	zsdmp = 1./(sdmptr * rday)
474	zbdmp = 1./(bdmptr * rday)
475	DO jn = 1, jptra
476	DO jk = 2, jpkm1
477	DO jj = 1, jpj
478	DO ji = 1, jpi
479	zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
480
481	! ... Decrease the value in the vicinity of the coast
482	restotr(ji,jj,jk,jn) = restotr(ji,jj,1,jn)*0.5 &
483	& * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )
484
485	! ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
486	restotr(ji,jj,jk,jn) = restotr(ji,jj,jk,jn) &
487	& * ( zbdmp + (zsdmp-zbdmp)*EXP(-fsdept(ji,jj,jk)/hdmptr) )
488	END DO
489	END DO
490	END DO
491	END DO
492
493	ENDIF
494
495
496	IF( cp_cfg == "orca" .AND. ( ndmptr > 0 .OR. ndmptr == -1 ) ) THEN
497
498	! ! =========================
499	! ! Med and Red Sea damping
500	! ! =========================
501	IF(lwp)WRITE(numout,*)
502	IF(lwp)WRITE(numout,*) ' ORCA configuration: Damping in Med and Red Seas'
503
504
505	zmrs(:,:) = 0.e0 ! damping term on the Med or Red Sea
506
507	SELECT CASE ( jp_cfg )
508	! ! =======================
509	CASE ( 4 ) ! ORCA_R4 configuration
510	! ! =======================
511
512	! Mediterranean Sea
513	ij0 = 50 ; ij1 = 56
514	ii0 = 81 ; ii1 = 91 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
515	ij0 = 50 ; ij1 = 55
516	ii0 = 75 ; ii1 = 80 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
517	ij0 = 52 ; ij1 = 53
518	ii0 = 70 ; ii1 = 74 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
519	! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
520	DO jk = 1, 17
521	zhfac (jk) = 0.5( 1.- COS( rpi(jk-1)/16. ) ) / rday
522	END DO
523	DO jk = 18, jpkm1
524	zhfac (jk) = 1./rday
525	END DO
526
527	! ! =======================
528	CASE ( 2 ) ! ORCA_R2 configuration
529	! ! =======================
530
531	! Mediterranean Sea
532	ij0 = 96 ; ij1 = 110
533	ii0 = 157 ; ii1 = 181 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
534	ij0 = 100 ; ij1 = 110
535	ii0 = 144 ; ii1 = 156 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
536	ij0 = 100 ; ij1 = 103
537	ii0 = 139 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
538	! Decrease before Gibraltar Strait
539	ij0 = 101 ; ij1 = 102
540	ii0 = 139 ; ii1 = 141 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
541	ii0 = 142 ; ii1 = 142 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
542	ii0 = 143 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
543	ii0 = 144 ; ii1 = 144 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0
544	! Red Sea
545	ij0 = 87 ; ij1 = 96
546	ii0 = 147 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
547	! Decrease before Bab el Mandeb Strait
548	ij0 = 91 ; ij1 = 91
549	ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0
550	ij0 = 90 ; ij1 = 90
551	ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
552	ij0 = 89 ; ij1 = 89
553	ii0 = 158 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
554	ij0 = 88 ; ij1 = 88
555	ii0 = 160 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
556	! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
557	DO jk = 1, 17
558	zhfac (jk) = 0.5( 1.- COS( rpi(jk-1)/16. ) ) / rday
559	END DO
560	DO jk = 18, jpkm1
561	zhfac (jk) = 1./rday
562	END DO
563
564	! ! =======================
565	CASE ( 05 ) ! ORCA_R05 configuration
566	! ! =======================
567
568	! Mediterranean Sea
569	ii0 = 568 ; ii1 = 574
570	ij0 = 324 ; ij1 = 333 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
571	ii0 = 575 ; ii1 = 658
572	ij0 = 314 ; ij1 = 366 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
573	! Black Sea (remaining part
574	ii0 = 641 ; ii1 = 651
575	ij0 = 367 ; ij1 = 372 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
576	! Decrease before Gibraltar Strait
577	ii0 = 324 ; ii1 = 333
578	ij0 = 565 ; ij1 = 565 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
579	ij0 = 566 ; ij1 = 566 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
580	ij0 = 567 ; ij1 = 567 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75
581	! Red Sea
582	ii0 = 641 ; ii1 = 665
583	ij0 = 270 ; ij1 = 310 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
584	! Decrease before Bab el Mandeb Strait
585	ii0 = 666 ; ii1 = 675
586	ij0 = 270 ; ij1 = 290
587	DO ji = mi0(ii0), mi1(ii1)
588	zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1 * ABS( FLOAT(ji - mi1(ii1)) )
589	END DO
590	zsdmp = 1./(sdmptr * rday)
591	zbdmp = 1./(bdmptr * rday)
592	DO jk = 1, jpk
593	zhfac (jk) = ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(1,1,jk)/hdmptr) )
594	END DO
595
596	! ! ========================
597	CASE ( 025 ) ! ORCA_R025 configuration
598	! ! ========================
599	IF(lwp) WRITE(numout,cform_err)
600	IF(lwp) WRITE(numout,*)' Not yet implemented in ORCA_R025'
601	nstop = nstop + 1
602
603	END SELECT
604
605	DO jn = 1, jptra
606	DO jk = 1, jpkm1
607	restotr(:,:,jk,jn) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * restotr(:,:,jk,jn)
608	END DO
609
610	! Mask resto array and set to 0 first and last levels
611	restotr(:,:, : ,jn) = restotr(:,:,:,jn) * tmask(:,:,:)
612	restotr(:,:, 1 ,jn) = 0.e0
613	restotr(:,:,jpk,jn) = 0.e0
614	END DO
615
616	ELSE
617	! ------------
618	! No damping
619	! ------------
620	IF(lwp) WRITE(numout,cform_err)
621	IF(lwp) WRITE(numout,*) 'Choose a correct value of ndmptr or DO NOT defined key_trcdmp'
622	nstop = nstop + 1
623	ENDIF
624
625	! ----------------------------
626	! Create Print damping array
627	! ----------------------------
628
629	! ndmpftr : = 1 create a damping.coeff NetCDF file
630
631	IF( ndmpftr == 1 ) THEN
632	DO jn = 1, jptra
633	IF(lwp) WRITE(numout,*) ' create damping.coeff.nc file ',jn
634	itime = 0
635	clname3 = 'damping.coeff'//ctrcnm(jn)
636	CALL ymds2ju( 0 , 1 , 1 , 0.e0 , zdate0 )
637	CALL restini( 'NONE', jpi , jpj , glamt, gphit, &
638	& jpk , gdept , clname3, itime, zdate0, &
639	& rdt , idmp )
640	CALL restput( idmp, 'Resto', jpi, jpj, jpk, 0 , restotr(:,:,:,jn) )
641	CALL restclo( idmp )
642	END DO
643	ENDIF
644
645
646	END SUBROUTINE trccof
647
648
649	SUBROUTINE cofdis ( pdct )
650	!!----------------------------------------------------------------------
651	!! * ROUTINE cofdis *
652	!!
653	!! ** Purpose : Compute the distance between ocean T-points and the
654	!! ocean model coastlines. Save the distance in a NetCDF file.
655	!!
656	!! ** Method : For each model level, the distance-to-coast is
657	!! computed as follows :
658	!! - The coastline is defined as the serie of U-,V-,F-points
659	!! that are at the ocean-land bound.
660	!! - For each ocean T-point, the distance-to-coast is then
661	!! computed as the smallest distance (on the sphere) between the
662	!! T-point and all the coastline points.
663	!! - For land T-points, the distance-to-coast is set to zero.
664	!! C A U T I O N : Computation not yet implemented in mpp case.
665	!!
666	!! ** Action : - pdct, distance to the coastline (argument)
667	!! - NetCDF file 'trc.dist.coast.nc'
668	!!
669	!! History :
670	!! 7.0 ! 01-02 (M. Imbard) Original code
671	!! 8.1 ! 01-02 (G. Madec, E. Durand)
672	!! 8.5 ! 02-08 (G. Madec, E. Durand) Free form, F90
673	!!----------------------------------------------------------------------
674	!! * Modules used
675	USE ioipsl
676
677	!! * Arguments
678	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
679	pdct ! distance to the coastline
680
681	!! * local declarations
682	INTEGER :: ji, jj, jk, jl ! dummy loop indices
683	INTEGER :: iju, ijt ! temporary integers
684	INTEGER :: icoast, itime
685	INTEGER :: &
686	icot ! logical unit for file distance to the coast
687	LOGICAL, DIMENSION(jpi,jpj) :: &
688	llcotu, llcotv, llcotf ! ???
689	CHARACTER (len=32) :: clname
690	REAL(wp) :: zdate0
691	REAL(wp), DIMENSION(jpi,jpj) :: &
692	zxt, zyt, zzt, & ! cartesian coordinates for T-points
693	zmask
694	REAL(wp), DIMENSION(3jpijpj) :: &
695	zxc, zyc, zzc, zdis ! temporary workspace
696	!!----------------------------------------------------------------------
697
698	! 0. Initialization
699	! -----------------
700	IF(lwp) WRITE(numout,*)
701	IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
702	IF(lwp) WRITE(numout,*) '~~~~~~'
703	IF(lwp) WRITE(numout,*)
704	IF( lk_mpp ) THEN
705	IF(lwp) WRITE(numout,cform_err)
706	IF(lwp) WRITE(numout,*) ' Computation not yet implemented with key_mpp_...'
707	IF(lwp) WRITE(numout,*) ' Rerun the code on another computer or '
708	IF(lwp) WRITE(numout,*) ' create the "dist.coast.nc" file using IDL'
709	nstop = nstop + 1
710	ENDIF
711
712	pdct(:,:,:) = 0.e0
713	zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )
714	zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )
715	zzt(:,:) = sin( rad * gphit(:,:) )
716
717
718	! 1. Loop on vertical levels
719	! --------------------------
720	! ! ===============
721	DO jk = 1, jpkm1 ! Horizontal slab
722	! ! ===============
723	! Define the coastline points (U, V and F)
724	DO jj = 2, jpjm1
725	DO ji = 2, jpim1
726	zmask(ji,jj) = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
727	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
728	llcotu(ji,jj) = ( tmask(ji,jj, jk) + tmask(ji+1,jj ,jk) == 1. )
729	llcotv(ji,jj) = ( tmask(ji,jj ,jk) + tmask(ji ,jj+1,jk) == 1. )
730	llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
731	END DO
732	END DO
733
734	! Lateral boundaries conditions
735	llcotu(:, 1 ) = umask(:, 2 ,jk) == 1
736	llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
737	llcotv(:, 1 ) = vmask(:, 2 ,jk) == 1
738	llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
739	llcotf(:, 1 ) = fmask(:, 2 ,jk) == 1
740	llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
741
742	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
743	llcotu( 1 ,:) = llcotu(jpim1,:)
744	llcotu(jpi,:) = llcotu( 2 ,:)
745	llcotv( 1 ,:) = llcotv(jpim1,:)
746	llcotv(jpi,:) = llcotv( 2 ,:)
747	llcotf( 1 ,:) = llcotf(jpim1,:)
748	llcotf(jpi,:) = llcotf( 2 ,:)
749	ELSE
750	llcotu( 1 ,:) = umask( 2 ,:,jk) == 1
751	llcotu(jpi,:) = umask(jpim1,:,jk) == 1
752	llcotv( 1 ,:) = vmask( 2 ,:,jk) == 1
753	llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
754	llcotf( 1 ,:) = fmask( 2 ,:,jk) == 1
755	llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
756	ENDIF
757	IF( nperio == 3 .OR. nperio == 4 ) THEN
758	DO ji = 1, jpim1
759	iju = jpi - ji + 1
760	llcotu(ji,jpj ) = llcotu(iju,jpj-2)
761	llcotf(ji,jpj-1) = llcotf(iju,jpj-2)
762	llcotf(ji,jpj ) = llcotf(iju,jpj-3)
763	END DO
764	DO ji = jpi/2, jpi-1
765	iju = jpi - ji + 1
766	llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
767	END DO
768	DO ji = 2, jpi
769	ijt = jpi - ji + 2
770	llcotv(ji,jpj-1) = llcotv(ijt,jpj-2)
771	llcotv(ji,jpj ) = llcotv(ijt,jpj-3)
772	END DO
773	ENDIF
774	IF( nperio == 5 .OR. nperio == 6 ) THEN
775	DO ji = 1, jpim1
776	iju = jpi - ji
777	llcotu(ji,jpj ) = llcotu(iju,jpj-1)
778	llcotf(ji,jpj ) = llcotf(iju,jpj-2)
779	END DO
780	DO ji = jpi/2, jpi-1
781	iju = jpi - ji
782	llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
783	END DO
784	DO ji = 1, jpi
785	ijt = jpi - ji + 1
786	llcotv(ji,jpj ) = llcotv(ijt,jpj-1)
787	END DO
788	DO ji = jpi/2+1, jpi
789	ijt = jpi - ji + 1
790	llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
791	END DO
792	ENDIF
793
794	! Compute cartesian coordinates of coastline points
795	! and the number of coastline points
796
797	icoast = 0
798	DO jj = 1, jpj
799	DO ji = 1, jpi
800	IF( llcotf(ji,jj) ) THEN
801	icoast = icoast + 1
802	zxc(icoast) = COS( radgphif(ji,jj) ) COS( rad*glamf(ji,jj) )
803	zyc(icoast) = COS( radgphif(ji,jj) ) SIN( rad*glamf(ji,jj) )
804	zzc(icoast) = SIN( rad*gphif(ji,jj) )
805	ENDIF
806	IF( llcotu(ji,jj) ) THEN
807	icoast = icoast+1
808	zxc(icoast) = COS( radgphiu(ji,jj) ) COS( rad*glamu(ji,jj) )
809	zyc(icoast) = COS( radgphiu(ji,jj) ) SIN( rad*glamu(ji,jj) )
810	zzc(icoast) = SIN( rad*gphiu(ji,jj) )
811	ENDIF
812	IF( llcotv(ji,jj) ) THEN
813	icoast = icoast+1
814	zxc(icoast) = COS( radgphiv(ji,jj) ) COS( rad*glamv(ji,jj) )
815	zyc(icoast) = COS( radgphiv(ji,jj) ) SIN( rad*glamv(ji,jj) )
816	zzc(icoast) = SIN( rad*gphiv(ji,jj) )
817	ENDIF
818	END DO
819	END DO
820
821	! Distance for the T-points
822
823	DO jj = 1, jpj
824	DO ji = 1, jpi
825	IF( tmask(ji,jj,jk) == 0. ) THEN
826	pdct(ji,jj,jk) = 0.
827	ELSE
828	DO jl = 1, icoast
829	zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2 &
830	+ ( zyt(ji,jj) - zyc(jl) )**2 &
831	+ ( zzt(ji,jj) - zzc(jl) )**2
832	END DO
833	pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
834	ENDIF
835	END DO
836	END DO
837	! ! ===============
838	END DO ! End of slab
839	! ! ===============
840
841
842	! 2. Create the distance to the coast file in NetCDF format
843	! ----------------------------------------------------------
844	clname = 'trc.dist.coast'
845	itime = 0
846	CALL ymds2ju( 0 , 1 , 1 , 0.e0 , zdate0 )
847	CALL restini( 'NONE', jpi , jpj , glamt, gphit , &
848	jpk , gdept , clname, itime, zdate0, &
849	rdt , icot )
850	CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
851	CALL restclo( icot )
852
853	END SUBROUTINE cofdis
854
855	#else
856	!!----------------------------------------------------------------------
857	!! Default key NO internal damping
858	!!----------------------------------------------------------------------
859	LOGICAL , PUBLIC, PARAMETER :: lk_trcdmp = .FALSE. !: internal damping flag
860	CONTAINS
861	SUBROUTINE trc_dmp( kt ) ! Empty routine
862	INTEGER, INTENT(in) :: kt
863	WRITE(,) 'trc_dmp: You should not have seen this print! error?', kt
864	END SUBROUTINE trc_dmp
865	#endif
866
867	!!======================================================================
868	END MODULE trcdmp

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

Download in other formats: