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.

trazdf_imp.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 15.6 KB

Rev	Line
[3]	1	MODULE trazdf_imp
[1438]	2	!!======================================================================
[457]	3	!! * MODULE trazdf_imp *
[2528]	4	!! Ocean tracers: vertical component of the tracer mixing trend
[1438]	5	!!======================================================================
	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 7.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction on tracer trend loops
	9	!! ! 1996-01 (G. Madec) statement function for e3
	10	!! ! 1997-05 (G. Madec) vertical component of isopycnal
	11	!! ! 1997-07 (G. Madec) geopotential diffusion in s-coord
	12	!! ! 2000-08 (G. Madec) double diffusive mixing
	13	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	14	!! 2.0 ! 2006-11 (G. Madec) New step reorganisation
	15	!! 3.2 ! 2009-03 (G. Madec) heat and salt content trends
[2528]	16	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC
[2602]	17	!! - ! 2011-02 (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition
[3]	18	!!----------------------------------------------------------------------
[1438]	19
	20	!!----------------------------------------------------------------------
[457]	21	!! tra_zdf_imp : Update the tracer trend with the diagonal vertical
	22	!! part of the mixing tensor.
[3]	23	!!----------------------------------------------------------------------
[457]	24	USE oce ! ocean dynamics and tracers variables
	25	USE dom_oce ! ocean space and time domain variables
	26	USE zdf_oce ! ocean vertical physics variables
[2528]	27	USE trc_oce ! share passive tracers/ocean variables
	28	USE domvvl ! variable volume
[216]	29	USE ldftra_oce ! ocean active tracers: lateral physics
[2528]	30	USE ldftra ! lateral mixing type
[457]	31	USE ldfslp ! lateral physics: slope of diffusion
	32	USE zdfddm ! ocean vertical physics: double diffusion
[2528]	33	USE traldf_iso_grif ! active tracers: Griffies operator
[3]	34	USE in_out_manager ! I/O manager
[457]	35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2715]	36	USE lib_mpp ! MPP library
[3]	37
	38	IMPLICIT NONE
	39	PRIVATE
	40
[1438]	41	PUBLIC tra_zdf_imp ! routine called by step.F90
[3]	42
[2602]	43	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
	44
[3211]	45	!! * Control permutation of array indices
	46	# include "oce_ftrans.h90"
	47	# include "dom_oce_ftrans.h90"
	48	# include "zdf_oce_ftrans.h90"
	49	# include "trc_oce_ftrans.h90"
	50	# include "domvvl_ftrans.h90"
	51	# include "ldftra_oce_ftrans.h90"
	52	# include "ldfslp_ftrans.h90"
	53	# include "zdfddm_ftrans.h90"
	54	# include "traldf_iso_grif_ftrans.h90"
	55
[3]	56	!! * Substitutions
	57	# include "domzgr_substitute.h90"
[457]	58	# include "ldftra_substitute.h90"
[3]	59	# include "zdfddm_substitute.h90"
[457]	60	# include "vectopt_loop_substitute.h90"
[3]	61	!!----------------------------------------------------------------------
[2528]	62	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]	63	!! $Id$
[2528]	64	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[457]	65	!!----------------------------------------------------------------------
[3]	66	CONTAINS
[2528]	67
	68	SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt )
[3]	69	!!----------------------------------------------------------------------
	70	!! * ROUTINE tra_zdf_imp *
	71	!!
[2602]	72	!! ** Purpose : Compute the after tracer through a implicit computation
	73	!! of the vertical tracer diffusion (including the vertical component
	74	!! of lateral mixing (only for 2nd order operator, for fourth order
	75	!! it is already computed and add to the general trend in traldf)
[3]	76	!!
[2602]	77	!! ** Method : The vertical diffusion of the tracer t is given by:
	78	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
[457]	79	!! It is computed using a backward time scheme (t=ta).
[2602]	80	!! If lk_zdfddm=T, use avs for salinity or for passive tracers
[3]	81	!! Surface and bottom boundary conditions: no diffusive flux on
	82	!! both tracers (bottom, applied through the masked field avt).
[2602]	83	!! If iso-neutral mixing, add to avt the contribution due to lateral mixing.
[3]	84	!!
[2602]	85	!! ** Action : - pta becomes the after tracer
[3]	86	!!---------------------------------------------------------------------
[2715]	87	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	88	USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace
	89	USE wrk_nemo, ONLY: zwi => wrk_3d_6 , zwt => wrk_3d_7 ! 3D workspace
[3211]	90
	91	!! DCSE_NEMO: Need additional directives for renamed module variables
	92	!FTRANS zwd zws :I :I :z
	93	!FTRANS zwi zwt :I :I :z
[2715]	94	!
[2528]	95	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	96	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	97	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	98	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	99
	100	!! DCSE_NEMO: This style defeats ftrans
	101	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	102	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	103
	104	!FTRANS ptb pta :I :I :z :
	105	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! before and now tracer fields
	106	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
[2715]	107	!
	108	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	109	REAL(wp) :: zrhs, ze3tb, ze3tn, ze3ta ! local scalars
[3]	110	!!---------------------------------------------------------------------
	111
[2715]	112	IF( wrk_in_use(3, 6,7) ) THEN
	113	CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.') ; RETURN
	114	ENDIF
	115
[2528]	116	IF( kt == nit000 ) THEN
[457]	117	IF(lwp)WRITE(numout,*)
[2528]	118	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
[457]	119	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
[2602]	120	!
	121	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
	122	ELSE ; r_vvl = 0._wp
	123	ENDIF
[457]	124	ENDIF
[2528]	125	!
[2602]	126	! ! ============= !
	127	DO jn = 1, kjpt ! tracer loop !
	128	! ! ============= !
[2528]	129	!
	130	! Matrix construction
[2602]	131	! --------------------
	132	! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
	133	!
[2715]	134	IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. &
[2602]	135	& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
	136	!
	137	! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
[3211]	138	#if defined key_z_first
	139	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
	142	zwt(ji,jj,1) = 0._wp
	143	DO jk = 2, jpk
	144	zwt(ji,jj,jk) = avt (ji,jj,jk)
	145	END DO
	146	END DO
	147	END DO
	148	ELSE
	149	DO jj = 1, jpj
	150	DO ji = 1, jpi
	151	zwt(ji,jj,1) = 0._wp
	152	DO jk = 2, jpk
	153	zwt(ji,jj,jk) = fsavs(ji,jj,jk)
	154	END DO
	155	END DO
	156	END DO
	157	ENDIF
	158	#else
[2602]	159	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt (:,:,2:jpk)
	160	ELSE ; zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
[2528]	161	ENDIF
[2602]	162	zwt(:,:,1) = 0._wp
[3211]	163	#endif
[2528]	164	!
	165	#if defined key_ldfslp
[2602]	166	! isoneutral diffusion: add the contribution
	167	IF( ln_traldf_grif ) THEN ! Griffies isoneutral diff
[3211]	168	#if defined key_z_first
	169	DO jj = 2, jpjm1
	170	DO ji = 2, jpim1
	171	DO jk = 2, jpkm1
	172	#else
[2528]	173	DO jk = 2, jpkm1
	174	DO jj = 2, jpjm1
	175	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	176	#endif
[2602]	177	zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
[2528]	178	END DO
	179	END DO
	180	END DO
[2602]	181	ELSE IF( l_traldf_rot ) THEN ! standard isoneutral diff
[3211]	182	#if defined key_z_first
	183	DO jj = 2, jpjm1
	184	DO ji = 2, jpim1
	185	DO jk = 2, jpkm1
	186	#else
[2528]	187	DO jk = 2, jpkm1
	188	DO jj = 2, jpjm1
	189	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	190	#endif
[2602]	191	zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk) &
	192	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	193	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
[2528]	194	END DO
	195	END DO
	196	END DO
	197	ENDIF
[592]	198	#endif
[2602]	199	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
[3211]	200	#if defined key_z_first
	201	DO jj = 2, jpjm1
	202	DO ji = 2, jpim1
	203	DO jk = 1, jpkm1
	204	ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	205	ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
	206	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	207	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	208	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	209	END DO
	210	#else
[2528]	211	DO jk = 1, jpkm1
	212	DO jj = 2, jpjm1
	213	DO ji = fs_2, fs_jpim1 ! vector opt.
[2602]	214	ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	215	ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
[2528]	216	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	217	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	218	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	219	END DO
	220	END DO
[3]	221	END DO
[3211]	222	#endif
[2528]	223	!
[2602]	224	!! Matrix inversion from the first level
	225	!!----------------------------------------------------------------------
	226	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	227	!
	228	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	229	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	230	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	231	! ( ... )( ... ) ( ... )
	232	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	233	!
	234	! m is decomposed in the product of an upper and lower triangular matrix.
	235	! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
	236	! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
	237	! and "superior" (above diagonal) components of the tridiagonal system.
	238	! The solution will be in the 4d array pta.
	239	! The 3d array zwt is used as a work space array.
	240	! En route to the solution pta is used a to evaluate the rhs and then
	241	! used as a work space array: its value is modified.
	242	!
[2528]	243	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
[2602]	244	! done once for all passive tracers (so included in the IF instruction)
[3211]	245	#if defined key_z_first
	246	zwt(ji,jj,1) = zwd(ji,jj,1)
	247	DO jk = 2, jpkm1
	248	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
	249	END DO
	250	END DO
	251	END DO
	252	#else
[2528]	253	DO jj = 2, jpjm1
	254	DO ji = fs_2, fs_jpim1
	255	zwt(ji,jj,1) = zwd(ji,jj,1)
	256	END DO
[457]	257	END DO
[2528]	258	DO jk = 2, jpkm1
	259	DO jj = 2, jpjm1
	260	DO ji = fs_2, fs_jpim1
	261	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
	262	END DO
	263	END DO
	264	END DO
[3211]	265	#endif
[2528]	266	!
	267	END IF
[2602]	268	!
[2528]	269	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
[3211]	270	#if defined key_z_first
[457]	271	DO jj = 2, jpjm1
[3211]	272	DO ji = 2, jpim1
	273	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
	274	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
	275	pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
	276	DO jk = 2, jpkm1
	277	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk)
	278	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk)
	279	zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side
	280	pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
	281	END DO
	282	#else
	283	DO jj = 2, jpjm1
[457]	284	DO ji = fs_2, fs_jpim1
[2602]	285	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
	286	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
[2528]	287	pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
[457]	288	END DO
	289	END DO
[2528]	290	DO jk = 2, jpkm1
	291	DO jj = 2, jpjm1
	292	DO ji = fs_2, fs_jpim1
[2602]	293	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk)
	294	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk)
[2528]	295	zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side
	296	pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
	297	END DO
[3]	298	END DO
	299	END DO
[3211]	300	#endif
[457]	301
[2602]	302	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
[3211]	303	#if defined key_z_first
	304	pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
	305	DO jk = jpk-2, 1, -1
	306	pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
	307	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	308	END DO
	309	END DO
	310	END DO
	311	#else
[457]	312	DO jj = 2, jpjm1
	313	DO ji = fs_2, fs_jpim1
[2528]	314	pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[3]	315	END DO
	316	END DO
[2528]	317	DO jk = jpk-2, 1, -1
	318	DO jj = 2, jpjm1
	319	DO ji = fs_2, fs_jpim1
	320	pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
	321	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	322	END DO
[457]	323	END DO
	324	END DO
[3211]	325	#endif
[2602]	326	! ! ================= !
	327	END DO ! end tracer loop !
	328	! ! ================= !
[1438]	329	!
[2715]	330	IF( wrk_not_released(3, 6,7) ) CALL ctl_stop('tra_zdf_imp: failed to release workspace arrays')
	331	!
[3]	332	END SUBROUTINE tra_zdf_imp
	333
	334	!!==============================================================================
	335	END MODULE trazdf_imp

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

Download in other formats: