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 branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2236

Last change on this file since 2236 was 2236, checked in by cetlod, 14 years ago
First guess of NEMO_v3.3
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 36.0 KB

Line
1	MODULE tradmp
2	!!======================================================================
3	!! * MODULE tradmp *
4	!! Ocean physics: internal restoring trend on active tracers (T and S)
5	!!======================================================================
6	!! History : OPA ! 1991-03 (O. Marti, G. Madec) Original code
7	!! ! 1992-06 (M. Imbard) doctor norme
8	!! ! 1996-01 (G. Madec) statement function for e3
9	!! ! 1997-05 (G. Madec) macro-tasked on jk-slab
10	!! ! 1998-07 (M. Imbard, G. Madec) ORCA version
11	!! 7.0 ! 2001-02 (M. Imbard) cofdis, Original code
12	!! 8.1 ! 2001-02 (G. Madec, E. Durand) cleaning
13	!! NEMO 1.0 ! 2002-08 (G. Madec, E. Durand) free form + modules
14	!! 3.2 ! 2009-08 (G. Madec, C. Talandier) DOCTOR norm for namelist parameter
15	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC
16	!!----------------------------------------------------------------------
17	#if defined key_tradmp \|\| defined key_esopa
18	!!----------------------------------------------------------------------
19	!! key_tradmp internal damping
20	!!----------------------------------------------------------------------
21	!! tra_dmp : update the tracer trend with the internal damping
22	!! tra_dmp_init : initialization, namlist read, parameters control
23	!! dtacof_zoom : restoring coefficient for zoom domain
24	!! dtacof : restoring coefficient for global domain
25	!! cofdis : compute the distance to the coastline
26	!!----------------------------------------------------------------------
27	USE oce ! ocean dynamics and tracers variables
28	USE dom_oce ! ocean space and time domain variables
29	USE trdmod_oce ! ocean space and time domain variables
30	USE trdtra ! ocean space and time domain variables
31	USE zdf_oce ! ocean vertical physics
32	USE phycst ! Define parameters for the routines
33	USE dtatem ! temperature data
34	USE dtasal ! salinity data
35	USE zdfmxl ! mixed layer depth
36	USE in_out_manager ! I/O manager
37	USE lib_mpp ! distribued memory computing
38	USE prtctl ! Print control
39
40	IMPLICIT NONE
41	PRIVATE
42
43	PUBLIC tra_dmp ! routine called by step.F90
44	PUBLIC tra_dmp_init ! routine called by opa.F90
45	PUBLIC dtacof ! routine called by tradmp.F90 and trcdmp.F90
46	PUBLIC dtacof_zoom ! routine called by tradmp.F90 and trcdmp.F90
47
48	#if ! defined key_agrif
49	LOGICAL, PUBLIC, PARAMETER :: lk_tradmp = .TRUE. !: internal damping flag
50	#else
51	LOGICAL, PUBLIC :: lk_tradmp = .TRUE. !: internal damping flag
52	#endif
53	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: strdmp !: damping salinity trend (psu/s)
54	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: ttrdmp !: damping temperature trend (Centigrade/s)
55	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: resto !: restoring coeff. on T and S (s-1)
56
57	! !!* Namelist namtra_dmp : T & S newtonian damping *
58	INTEGER :: nn_hdmp = -1 ! = 0/-1/'latitude' for damping over T and S
59	INTEGER :: nn_zdmp = 0 ! = 0/1/2 flag for damping in the mixed layer
60	REAL(wp) :: rn_surf = 50. ! surface time scale for internal damping [days]
61	REAL(wp) :: rn_bot = 360. ! bottom time scale for internal damping [days]
62	REAL(wp) :: rn_dep = 800. ! depth of transition between rn_surf and rn_bot [meters]
63	INTEGER :: nn_file = 2 ! = 1 create a damping.coeff NetCDF file
64
65	!! * Substitutions
66	# include "domzgr_substitute.h90"
67	# include "vectopt_loop_substitute.h90"
68	!!----------------------------------------------------------------------
69	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
70	!! $Id$
71	!! Software governed by the CeCILL licence (NEMOGCM/License_CeCILL.txt)
72	!!----------------------------------------------------------------------
73
74	CONTAINS
75
76	SUBROUTINE tra_dmp( kt )
77	!!----------------------------------------------------------------------
78	!! * ROUTINE tra_dmp *
79	!!
80	!! ** Purpose : Compute the tracer trend due to a newtonian damping
81	!! of the tracer field towards given data field and add it to the
82	!! general tracer trends.
83	!!
84	!! ** Method : Newtonian damping towards t_dta and s_dta computed
85	!! and add to the general tracer trends:
86	!! ta = ta + resto * (t_dta - tb)
87	!! sa = sa + resto * (s_dta - sb)
88	!! The trend is computed either throughout the water column
89	!! (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
90	!! below the well mixed layer (nlmdmp=2)
91	!!
92	!! ** Action : - (ta,sa) tracer trends updated with the damping trend
93	!!----------------------------------------------------------------------
94	INTEGER, INTENT(in) :: kt ! ocean time-step index
95	!!
96	REAL(wp) :: zta, zsa ! temporary scalars
97	INTEGER :: ji, jj, jk ! dummy loop indices
98	REAL(wp) :: zta, zsa
99	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
100	!!----------------------------------------------------------------------
101
102	IF( l_trdtra ) THEN !* Save ta and sa trends
103	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
104	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
105	ENDIF
106
107	SELECT CASE ( nn_zdmp )
108	!
109	CASE( 0 ) !== newtonian damping throughout the water column ==!
110	DO jk = 1, jpkm1
111	DO jj = 2, jpjm1
112	DO ji = fs_2, fs_jpim1 ! vector opt.
113	zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
114	zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
115	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
116	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
117	! save the salinity trend (used in asmtrj)
118	strdmp(ji,jj,jk) = zsa
119	ttrdmp(ji,jj,jk) = zta
120	END DO
121	END DO
122	END DO
123	!
124	CASE ( 1 ) !== no damping in the turbocline (avt > 5 cm2/s) ==!
125	DO jk = 1, jpkm1
126	DO jj = 2, jpjm1
127	DO ji = fs_2, fs_jpim1 ! vector opt.
128	IF( avt(ji,jj,jk) <= 5.e-4 ) THEN
129	zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
130	zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
131	ELSE
132	zta = 0.e0
133	zsa = 0.e0
134	ENDIF
135	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
136	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
137	! save the salinity trend (used in asmtrj)
138	strdmp(ji,jj,jk) = zsa
139	ttrdmp(ji,jj,jk) = zta
140	END DO
141	END DO
142	END DO
143	!
144	CASE ( 2 ) !== no damping in the mixed layer ==!
145	DO jk = 1, jpkm1
146	DO jj = 2, jpjm1
147	DO ji = fs_2, fs_jpim1 ! vector opt.
148	IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
149	zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
150	zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
151	ELSE
152	zta = 0.e0
153	zsa = 0.e0
154	ENDIF
155	! add the trends to the general tracer trends
156	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
157	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
158	! save the salinity trend (used in asmtrj)
159	strdmp(ji,jj,jk) = zsa
160	ttrdmp(ji,jj,jk) = zta
161	END DO
162	END DO
163	END DO
164	!
165	END SELECT
166
167	IF( l_trdtra ) THEN ! trend diagnostic
168	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
169	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
170	CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ztrdt )
171	CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, ztrds )
172	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
173	ENDIF
174	! ! Control print
175	IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' dmp - Ta: ', mask1=tmask, &
176	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
177	!
178	END SUBROUTINE tra_dmp
179
180
181	SUBROUTINE tra_dmp_init
182	!!----------------------------------------------------------------------
183	!! * ROUTINE tra_dmp_init *
184	!!
185	!! ** Purpose : Initialization for the newtonian damping
186	!!
187	!! ** Method : read the nammbf namelist and check the parameters
188	!!----------------------------------------------------------------------
189	NAMELIST/namtra_dmp/ nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file
190	!!----------------------------------------------------------------------
191
192	REWIND ( numnam ) ! Read Namelist namtra_dmp : temperature and salinity damping term
193	READ ( numnam, namtra_dmp )
194	IF( lzoom ) nn_zdmp = 0 ! restoring to climatology at closed north or south boundaries
195
196	IF(lwp) THEN ! Namelist print
197	WRITE(numout,*)
198	WRITE(numout,*) 'tra_dmp : T and S newtonian damping'
199	WRITE(numout,*) '~~~~~~~'
200	WRITE(numout,*) ' Namelist namtra_dmp : set damping parameter'
201	WRITE(numout,*) ' T and S damping option nn_hdmp = ', nn_hdmp
202	WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp, '(zoom: forced to 0)'
203	WRITE(numout,*) ' surface time scale (days) rn_surf = ', rn_surf
204	WRITE(numout,*) ' bottom time scale (days) rn_bot = ', rn_bot
205	WRITE(numout,*) ' depth of transition (meters) rn_dep = ', rn_dep
206	WRITE(numout,*) ' create a damping.coeff file nn_file = ', nn_file
207	ENDIF
208
209	SELECT CASE ( nn_hdmp )
210	CASE ( -1 ) ; IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only'
211	CASE ( 1:90 ) ; IF(lwp) WRITE(numout,*) ' tracer damping poleward of', nn_hdmp, ' degrees'
212	CASE DEFAULT
213	WRITE(ctmp1,*) ' bad flag value for nn_hdmp = ', nn_hdmp
214	CALL ctl_stop(ctmp1)
215	END SELECT
216
217	SELECT CASE ( nn_zdmp )
218	CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column'
219	CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline (avt > 5 cm2/s)'
220	CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer'
221	CASE DEFAULT
222	WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp
223	CALL ctl_stop(ctmp1)
224	END SELECT
225
226	IF( .NOT.lk_dtasal .OR. .NOT.lk_dtatem ) &
227	& CALL ctl_stop( 'no temperature and/or salinity data define key_dtatem and key_dtasal' )
228
229	strdmp(:,:,:) = 0.e0 ! internal damping salinity trend (used in asmtrj)
230	ttrdmp(:,:,:) = 0.e0
231	! ! Damping coefficients initialization
232	IF( lzoom ) THEN ; CALL dtacof_zoom( resto )
233	ELSE ; CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, &
234	& nn_file, 'TRA' , resto )
235	ENDIF
236	!
237	END SUBROUTINE tra_dmp_init
238
239
240	SUBROUTINE dtacof_zoom( presto )
241	!!----------------------------------------------------------------------
242	!! * ROUTINE dtacof_zoom *
243	!!
244	!! ** Purpose : Compute the damping coefficient for zoom domain
245	!!
246	!! ** Method : - set along closed boundary due to zoom a damping over
247	!! 6 points with a max time scale of 5 days.
248	!! - ORCA arctic/antarctic zoom: set the damping along
249	!! south/north boundary over a latitude strip.
250	!!
251	!! ** Action : - resto, the damping coeff. for T and S
252	!!----------------------------------------------------------------------
253	!! * Arguments
254	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto !: restoring coeff. (s-1)
255	!
256	INTEGER :: ji, jj, jk, jn ! dummy loop indices
257	REAL(wp) :: zlat, zlat0, zlat1, zlat2 ! temporary scalar
258	REAL(wp), DIMENSION(6) :: zfact ! temporary workspace
259	!!----------------------------------------------------------------------
260
261	zfact(1) = 1.
262	zfact(2) = 1.
263	zfact(3) = 11./12.
264	zfact(4) = 8./12.
265	zfact(5) = 4./12.
266	zfact(6) = 1./12.
267	zfact(:) = zfact(:) / ( 5. * rday ) ! 5 days max restoring time scale
268
269	presto(:,:,:) = 0.e0
270
271	! damping along the forced closed boundary over 6 grid-points
272	DO jn = 1, 6
273	IF( lzoom_w ) presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : ) = zfact(jn) ! west closed
274	IF( lzoom_s ) presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : ) = zfact(jn) ! south closed
275	IF( lzoom_e ) presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn) ! east closed
276	IF( lzoom_n ) presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn) ! north closed
277	END DO
278
279	! ! ====================================================
280	IF( lzoom_arct .AND. lzoom_anta ) THEN ! ORCA configuration : arctic zoom or antarctic zoom
281	! ! ====================================================
282	IF(lwp) WRITE(numout,*)
283	IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' dtacof_zoom : ORCA Arctic zoom'
284	IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' dtacof_zoom : ORCA Antarctic zoom'
285	IF(lwp) WRITE(numout,*)
286	!
287	! ! Initialization :
288	presto(:,:,:) = 0.e0
289	zlat0 = 10. ! zlat0 : latitude strip where resto decreases
290	zlat1 = 30. ! zlat1 : resto = 1 before zlat1
291	zlat2 = zlat1 + zlat0 ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
292
293	DO jk = 2, jpkm1 ! Compute arrays resto ; value for internal damping : 5 days
294	DO jj = 1, jpj
295	DO ji = 1, jpi
296	zlat = ABS( gphit(ji,jj) )
297	IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
298	presto(ji,jj,jk) = 0.5 * ( 1./(5.rday) ) ( 1. - cos(rpi*(zlat2-zlat)/zlat0) )
299	ELSEIF( zlat < zlat1 ) THEN
300	presto(ji,jj,jk) = 1./(5.*rday)
301	ENDIF
302	END DO
303	END DO
304	END DO
305	!
306	ENDIF
307	! ! Mask resto array
308	presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)
309	!
310	END SUBROUTINE dtacof_zoom
311
312
313	SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep, &
314	& kn_file, cdtype , presto )
315	!!----------------------------------------------------------------------
316	!! * ROUTINE dtacof *
317	!!
318	!! ** Purpose : Compute the damping coefficient
319	!!
320	!! ** Method : Arrays defining the damping are computed for each grid
321	!! point for temperature and salinity (resto)
322	!! Damping depends on distance to coast, depth and latitude
323	!!
324	!! ** Action : - resto, the damping coeff. for T and S
325	!!----------------------------------------------------------------------
326	USE iom
327	USE ioipsl
328	!! * Arguments
329	INTEGER , INTENT(in ) :: kn_hdmp !: damping option
330	REAL(wp) , INTENT(in ) :: pn_surf !: surface time scale (days)
331	REAL(wp) , INTENT(in ) :: pn_bot !: bottom time scale (days)
332	REAL(wp) , INTENT(in ) :: pn_dep !: depth of transition (meters)
333	INTEGER , INTENT(in ) :: kn_file !: save the damping coef on a file or not
334	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
335	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: presto !: restoring coeff. (s-1)
336	!
337	INTEGER :: ji, jj, jk ! dummy loop indices
338	INTEGER :: ii0, ii1, ij0, ij1 ! - -
339	INTEGER :: inum0 ! logical unit for file restoring damping term
340	INTEGER :: icot ! logical unit for file distance to the coast
341	REAL(wp) :: zinfl, zlon ! temporary scalars
342	REAL(wp) :: zlat, zlat0, zlat1, zlat2 ! - -
343	REAL(wp) :: zsdmp, zbdmp ! - -
344	REAL(wp), DIMENSION(jpk) :: zhfac
345	REAL(wp), DIMENSION(jpi,jpj) :: zmrs
346	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdct
347	CHARACTER(len=20) :: cfile
348	!!----------------------------------------------------------------------
349
350	! ====================================
351	! ORCA configuration : global domain
352	! ====================================
353
354	IF(lwp) WRITE(numout,*)
355	IF(lwp) WRITE(numout,*) ' dtacof : Global domain of ORCA'
356	IF(lwp) WRITE(numout,*) ' ------------------------------'
357
358	! ... Initialization :
359	presto(:,:,:) = 0.e0
360	!
361	IF( kn_hdmp > 0 ) THEN ! Damping poleward of 'nn_hdmp' degrees !
362	! !-----------------------------------------!
363	IF(lwp) WRITE(numout,*)
364	IF(lwp) WRITE(numout,*) ' Damping poleward of ', kn_hdmp,' deg.'
365	!
366	CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )
367	!
368	IF( icot > 0 ) THEN ! distance-to-coast read in file
369	CALL iom_get ( icot, jpdom_data, 'Tcoast', zdct )
370	CALL iom_close( icot )
371	ELSE ! distance-to-coast computed and saved in file (output in zdct)
372	CALL cofdis( zdct )
373	ENDIF
374
375	! ! Compute arrays resto
376	zinfl = 1000.e3 ! distance of influence for damping term
377	zlat0 = 10. ! latitude strip where resto decreases
378	zlat1 = REAL( kn_hdmp ) ! resto = 0 between -zlat1 and zlat1
379	zlat2 = zlat1 + zlat0 ! resto increases from 0 to 1 between \|zlat1\| and \|zlat2\|
380
381	DO jj = 1, jpj
382	DO ji = 1, jpi
383	zlat = ABS( gphit(ji,jj) )
384	IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
385	presto(ji,jj,1) = 0.5 * ( 1. - cos(rpi*(zlat-zlat1)/zlat0 ) )
386	ELSEIF ( zlat > zlat2 ) THEN
387	presto(ji,jj,1) = 1.
388	ENDIF
389	END DO
390	END DO
391
392	IF ( kn_hdmp == 20 ) THEN ! North Indian ocean (20N/30N x 45E/100E) : resto=0
393	DO jj = 1, jpj
394	DO ji = 1, jpi
395	zlat = gphit(ji,jj)
396	zlon = MOD( glamt(ji,jj), 360. )
397	IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45. < zlon .AND. zlon < 100. ) THEN
398	presto(ji,jj,1) = 0.e0
399	ENDIF
400	END DO
401	END DO
402	ENDIF
403
404	zsdmp = 1./(pn_surf * rday)
405	zbdmp = 1./(pn_bot * rday)
406	DO jk = 2, jpkm1
407	DO jj = 1, jpj
408	DO ji = 1, jpi
409	zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
410	! ... Decrease the value in the vicinity of the coast
411	presto(ji,jj,jk) = presto(ji,jj,1) * 0.5 * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )
412	! ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
413	presto(ji,jj,jk) = presto(ji,jj,jk) * ( zbdmp + (zsdmp-zbdmp)*EXP(-fsdept(ji,jj,jk)/pn_dep) )
414	END DO
415	END DO
416	END DO
417	!
418	ENDIF
419
420
421	IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN
422
423	! ! =========================
424	! ! Med and Red Sea damping
425	! ! =========================
426	IF(lwp)WRITE(numout,*)
427	IF(lwp)WRITE(numout,*) ' ORCA configuration: Damping in Med and Red Seas'
428
429
430	zmrs(:,:) = 0.e0 ! damping term on the Med or Red Sea
431
432	SELECT CASE ( jp_cfg )
433	! ! =======================
434	CASE ( 4 ) ! ORCA_R4 configuration
435	! ! =======================
436	! Mediterranean Sea
437	ij0 = 50 ; ij1 = 56
438	ii0 = 81 ; ii1 = 91 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
439	ij0 = 50 ; ij1 = 55
440	ii0 = 75 ; ii1 = 80 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
441	ij0 = 52 ; ij1 = 53
442	ii0 = 70 ; ii1 = 74 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
443	! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
444	DO jk = 1, 17
445	zhfac (jk) = 0.5( 1.- COS( rpi(jk-1)/16. ) ) / rday
446	END DO
447	DO jk = 18, jpkm1
448	zhfac (jk) = 1./rday
449	END DO
450	! ! =======================
451	CASE ( 2 ) ! ORCA_R2 configuration
452	! ! =======================
453	! Mediterranean Sea
454	ij0 = 96 ; ij1 = 110
455	ii0 = 157 ; ii1 = 181 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
456	ij0 = 100 ; ij1 = 110
457	ii0 = 144 ; ii1 = 156 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
458	ij0 = 100 ; ij1 = 103
459	ii0 = 139 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
460	! Decrease before Gibraltar Strait
461	ij0 = 101 ; ij1 = 102
462	ii0 = 139 ; ii1 = 141 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
463	ii0 = 142 ; ii1 = 142 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
464	ii0 = 143 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
465	ii0 = 144 ; ii1 = 144 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0
466	! Red Sea
467	ij0 = 87 ; ij1 = 96
468	ii0 = 147 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
469	! Decrease before Bab el Mandeb Strait
470	ij0 = 91 ; ij1 = 91
471	ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0
472	ij0 = 90 ; ij1 = 90
473	ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
474	ij0 = 89 ; ij1 = 89
475	ii0 = 158 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
476	ij0 = 88 ; ij1 = 88
477	ii0 = 160 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
478	! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
479	DO jk = 1, 17
480	zhfac (jk) = 0.5( 1.- COS( rpi(jk-1)/16. ) ) / rday
481	END DO
482	DO jk = 18, jpkm1
483	zhfac (jk) = 1./rday
484	END DO
485	! ! =======================
486	CASE ( 05 ) ! ORCA_R05 configuration
487	! ! =======================
488	! Mediterranean Sea
489	ii0 = 568 ; ii1 = 574
490	ij0 = 324 ; ij1 = 333 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
491	ii0 = 575 ; ii1 = 658
492	ij0 = 314 ; ij1 = 366 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
493	! Black Sea (remaining part
494	ii0 = 641 ; ii1 = 651
495	ij0 = 367 ; ij1 = 372 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
496	! Decrease before Gibraltar Strait
497	ij0 = 324 ; ij1 = 333
498	ii0 = 565 ; ii1 = 565 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
499	ii0 = 566 ; ii1 = 566 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
500	ii0 = 567 ; ii1 = 567 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75
501	! Red Sea
502	ii0 = 641 ; ii1 = 665
503	ij0 = 270 ; ij1 = 310 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
504	! Decrease before Bab el Mandeb Strait
505	ii0 = 666 ; ii1 = 675
506	ij0 = 270 ; ij1 = 290
507	DO ji = mi0(ii0), mi1(ii1)
508	zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1 * ABS( FLOAT(ji - mi1(ii1)) )
509	END DO
510	zsdmp = 1./(pn_surf * rday)
511	zbdmp = 1./(pn_bot * rday)
512	DO jk = 1, jpk
513	zhfac (jk) = ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(1,1,jk)/pn_dep) )
514	END DO
515	! ! ========================
516	CASE ( 025 ) ! ORCA_R025 configuration
517	! ! ========================
518	CALL ctl_stop( ' Not yet implemented in ORCA_R025' )
519	!
520	END SELECT
521
522	DO jk = 1, jpkm1
523	presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * presto(:,:,jk)
524	END DO
525
526	! Mask resto array and set to 0 first and last levels
527	presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
528	presto(:,:, 1 ) = 0.e0
529	presto(:,:,jpk) = 0.e0
530	! !--------------------!
531	ELSE ! No damping !
532	! !--------------------!
533	CALL ctl_stop( 'Choose a correct value of nn_hdmp or DO NOT defined key_tradmp' )
534	ENDIF
535
536	! !--------------------------------!
537	IF( kn_file == 1 ) THEN ! save damping coef. in a file !
538	! !--------------------------------!
539	IF(lwp) WRITE(numout,*) ' create damping.coeff.nc file'
540	IF( cdtype == 'TRA' ) cfile = 'damping.coeff'
541	IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'
542	cfile = TRIM( cfile )
543	CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
544	CALL iom_rstput( 0, 0, inum0, 'Resto', presto )
545	CALL iom_close ( inum0 )
546	ENDIF
547	!
548	END SUBROUTINE dtacof
549
550
551	SUBROUTINE cofdis( pdct )
552	!!----------------------------------------------------------------------
553	!! * ROUTINE cofdis *
554	!!
555	!! ** Purpose : Compute the distance between ocean T-points and the
556	!! ocean model coastlines. Save the distance in a NetCDF file.
557	!!
558	!! ** Method : For each model level, the distance-to-coast is
559	!! computed as follows :
560	!! - The coastline is defined as the serie of U-,V-,F-points
561	!! that are at the ocean-land bound.
562	!! - For each ocean T-point, the distance-to-coast is then
563	!! computed as the smallest distance (on the sphere) between the
564	!! T-point and all the coastline points.
565	!! - For land T-points, the distance-to-coast is set to zero.
566	!! C A U T I O N : Computation not yet implemented in mpp case.
567	!!
568	!! ** Action : - pdct, distance to the coastline (argument)
569	!! - NetCDF file 'dist.coast.nc'
570	!!----------------------------------------------------------------------
571	USE ioipsl ! IOipsl librairy
572	!!
573	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: pdct ! distance to the coastline
574	!!
575	INTEGER :: ji, jj, jk, jl ! dummy loop indices
576	INTEGER :: iju, ijt ! temporary integers
577	INTEGER :: icoast, itime
578	INTEGER :: icot ! logical unit for file distance to the coast
579	LOGICAL, DIMENSION(jpi,jpj) :: llcotu, llcotv, llcotf ! ???
580	CHARACTER (len=32) :: clname
581	REAL(wp) :: zdate0
582	REAL(wp), DIMENSION(jpi,jpj) :: zxt, zyt, zzt, zmask ! cartesian coordinates for T-points
583	REAL(wp), DIMENSION(3jpijpj) :: zxc, zyc, zzc, zdis ! temporary workspace
584	!!----------------------------------------------------------------------
585
586	! 0. Initialization
587	! -----------------
588	IF(lwp) WRITE(numout,*)
589	IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
590	IF(lwp) WRITE(numout,*) '~~~~~~'
591	IF(lwp) WRITE(numout,*)
592	IF( lk_mpp ) &
593	& CALL ctl_stop(' Computation not yet implemented with key_mpp_...', &
594	& ' Rerun the code on another computer or ', &
595	& ' create the "dist.coast.nc" file using IDL' )
596
597	pdct(:,:,:) = 0.e0
598	zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )
599	zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )
600	zzt(:,:) = sin( rad * gphit(:,:) )
601
602
603	! 1. Loop on vertical levels
604	! --------------------------
605	! ! ===============
606	DO jk = 1, jpkm1 ! Horizontal slab
607	! ! ===============
608	! Define the coastline points (U, V and F)
609	DO jj = 2, jpjm1
610	DO ji = 2, jpim1
611	zmask(ji,jj) = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
612	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
613	llcotu(ji,jj) = ( tmask(ji,jj, jk) + tmask(ji+1,jj ,jk) == 1. )
614	llcotv(ji,jj) = ( tmask(ji,jj ,jk) + tmask(ji ,jj+1,jk) == 1. )
615	llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
616	END DO
617	END DO
618
619	! Lateral boundaries conditions
620	llcotu(:, 1 ) = umask(:, 2 ,jk) == 1
621	llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
622	llcotv(:, 1 ) = vmask(:, 2 ,jk) == 1
623	llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
624	llcotf(:, 1 ) = fmask(:, 2 ,jk) == 1
625	llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
626
627	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
628	llcotu( 1 ,:) = llcotu(jpim1,:)
629	llcotu(jpi,:) = llcotu( 2 ,:)
630	llcotv( 1 ,:) = llcotv(jpim1,:)
631	llcotv(jpi,:) = llcotv( 2 ,:)
632	llcotf( 1 ,:) = llcotf(jpim1,:)
633	llcotf(jpi,:) = llcotf( 2 ,:)
634	ELSE
635	llcotu( 1 ,:) = umask( 2 ,:,jk) == 1
636	llcotu(jpi,:) = umask(jpim1,:,jk) == 1
637	llcotv( 1 ,:) = vmask( 2 ,:,jk) == 1
638	llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
639	llcotf( 1 ,:) = fmask( 2 ,:,jk) == 1
640	llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
641	ENDIF
642	IF( nperio == 3 .OR. nperio == 4 ) THEN
643	DO ji = 1, jpim1
644	iju = jpi - ji + 1
645	llcotu(ji,jpj ) = llcotu(iju,jpj-2)
646	llcotf(ji,jpjm1) = llcotf(iju,jpj-2)
647	llcotf(ji,jpj ) = llcotf(iju,jpj-3)
648	END DO
649	DO ji = jpi/2, jpim1
650	iju = jpi - ji + 1
651	llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
652	END DO
653	DO ji = 2, jpi
654	ijt = jpi - ji + 2
655	llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)
656	llcotv(ji,jpj ) = llcotv(ijt,jpj-3)
657	END DO
658	ENDIF
659	IF( nperio == 5 .OR. nperio == 6 ) THEN
660	DO ji = 1, jpim1
661	iju = jpi - ji
662	llcotu(ji,jpj ) = llcotu(iju,jpjm1)
663	llcotf(ji,jpj ) = llcotf(iju,jpj-2)
664	END DO
665	DO ji = jpi/2, jpim1
666	iju = jpi - ji
667	llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
668	END DO
669	DO ji = 1, jpi
670	ijt = jpi - ji + 1
671	llcotv(ji,jpj ) = llcotv(ijt,jpjm1)
672	END DO
673	DO ji = jpi/2+1, jpi
674	ijt = jpi - ji + 1
675	llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
676	END DO
677	ENDIF
678
679	! Compute cartesian coordinates of coastline points
680	! and the number of coastline points
681
682	icoast = 0
683	DO jj = 1, jpj
684	DO ji = 1, jpi
685	IF( llcotf(ji,jj) ) THEN
686	icoast = icoast + 1
687	zxc(icoast) = COS( radgphif(ji,jj) ) COS( rad*glamf(ji,jj) )
688	zyc(icoast) = COS( radgphif(ji,jj) ) SIN( rad*glamf(ji,jj) )
689	zzc(icoast) = SIN( rad*gphif(ji,jj) )
690	ENDIF
691	IF( llcotu(ji,jj) ) THEN
692	icoast = icoast+1
693	zxc(icoast) = COS( radgphiu(ji,jj) ) COS( rad*glamu(ji,jj) )
694	zyc(icoast) = COS( radgphiu(ji,jj) ) SIN( rad*glamu(ji,jj) )
695	zzc(icoast) = SIN( rad*gphiu(ji,jj) )
696	ENDIF
697	IF( llcotv(ji,jj) ) THEN
698	icoast = icoast+1
699	zxc(icoast) = COS( radgphiv(ji,jj) ) COS( rad*glamv(ji,jj) )
700	zyc(icoast) = COS( radgphiv(ji,jj) ) SIN( rad*glamv(ji,jj) )
701	zzc(icoast) = SIN( rad*gphiv(ji,jj) )
702	ENDIF
703	END DO
704	END DO
705
706	! Distance for the T-points
707
708	DO jj = 1, jpj
709	DO ji = 1, jpi
710	IF( tmask(ji,jj,jk) == 0. ) THEN
711	pdct(ji,jj,jk) = 0.
712	ELSE
713	DO jl = 1, icoast
714	zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2 &
715	& + ( zyt(ji,jj) - zyc(jl) )**2 &
716	& + ( zzt(ji,jj) - zzc(jl) )**2
717	END DO
718	pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
719	ENDIF
720	END DO
721	END DO
722	! ! ===============
723	END DO ! End of slab
724	! ! ===============
725
726
727	! 2. Create the distance to the coast file in NetCDF format
728	! ----------------------------------------------------------
729	clname = 'dist.coast'
730	itime = 0
731	CALL ymds2ju( 0 , 1 , 1 , 0.e0 , zdate0 )
732	CALL restini( 'NONE', jpi , jpj , glamt, gphit , &
733	& jpk , gdept_0, clname, itime, zdate0, &
734	& rdt , icot )
735	CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
736	CALL restclo( icot )
737
738	END SUBROUTINE cofdis
739
740	#else
741	!!----------------------------------------------------------------------
742	!! Default key NO internal damping
743	!!----------------------------------------------------------------------
744	LOGICAL , PUBLIC, PARAMETER :: lk_tradmp = .FALSE. !: internal damping flag
745	CONTAINS
746	SUBROUTINE tra_dmp( kt ) ! Empty routine
747	WRITE(,) 'tra_dmp: You should not have seen this print! error?', kt
748	END SUBROUTINE tra_dmp
749	SUBROUTINE tra_dmp_init ! Empty routine
750	END SUBROUTINE tra_dmp_init
751	#endif
752
753	!!======================================================================
754	END MODULE tradmp

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

Download in other formats: