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 @ 1076

Last change on this file since 1076 was 1076, checked in by cetlod, 16 years ago
update trcdmp.F90 module, see ticket:190
Property svn:executable set to ``*
File size: 34.2 KB

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

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

Download in other formats: