Context Navigation

tradmp.F90 @ 5098

Last change on this file since 5098 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

Property svn:keywords set to Id

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

Download in other formats: