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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2330

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

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

Download in other formats: