Context Navigation

tradmp.F90 @ 2416

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

Download in other formats: