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 @ 2787

Last change on this file since 2787 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 37.7 KB

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

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

Download in other formats: