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

Last change on this file since 1838 was 1601, checked in by ctlod, 15 years ago
Doctor naming of OPA namelist variables , see ticket: #526
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 33.0 KB

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

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

Download in other formats: