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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

tradmp.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 3787

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

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

Download in other formats: