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

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

tradmp.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 4328

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

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

Download in other formats: