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

Last change on this file since 941 was 941, checked in by cetlod, 16 years ago
phasing the passive tracer transport module to the new version of NEMO, see ticket 143
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 dta_trc( 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: