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.

tradmp.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/tradmp.F90 @ 225

Last change on this file since 225 was 225, checked in by opalod, 19 years ago
CT : UPDATE169 : correct indices location for damping coeficient close to the Gibraltar strait in the ORCA 0.5� configuration
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 37.1 KB

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

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

Download in other formats: