Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

tradmp.F90 @ 2578

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/TAM_V3_2_2/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2578

Download in other formats: