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/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2789

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

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

Download in other formats: