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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 6043

Last change on this file since 6043 was 6043, checked in by timgraham, 8 years ago
Merged head of trunk into branch
Property svn:keywords set to `Id`
File size: 11.7 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	!!----------------------------------------------------------------------
[6043]	21	!! tra_zdf_imp : Update the tracer trend with the diagonal vertical part of the mixing tensor.
[3]	22	!!----------------------------------------------------------------------
[6043]	23	USE oce ! ocean dynamics and tracers variables
	24	USE dom_oce ! ocean space and time domain variables
	25	USE zdf_oce ! ocean vertical physics variables
	26	USE trc_oce ! share passive tracers/ocean variables
	27	USE domvvl ! variable volume
	28	USE ldftra ! lateral mixing type
	29	USE ldfslp ! lateral physics: slope of diffusion
	30	USE zdfddm ! ocean vertical physics: double diffusion
	31	USE traldf_triad ! active tracers: Method of Stabilizing Correction
	32	!
	33	USE in_out_manager ! I/O manager
	34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	35	USE lib_mpp ! MPP library
	36	USE wrk_nemo ! Memory Allocation
	37	USE timing ! Timing
[3]	38
	39	IMPLICIT NONE
	40	PRIVATE
	41
[1438]	42	PUBLIC tra_zdf_imp ! routine called by step.F90
[3]	43
[2602]	44	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
	45
[3]	46	!! * Substitutions
	47	# include "domzgr_substitute.h90"
	48	# include "zdfddm_substitute.h90"
[457]	49	# include "vectopt_loop_substitute.h90"
[3]	50	!!----------------------------------------------------------------------
[6043]	51	!! NEMO/OPA 3.7 , NEMO Consortium (2015)
[1156]	52	!! $Id$
[2528]	53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[457]	54	!!----------------------------------------------------------------------
[3]	55	CONTAINS
[2528]	56
[3294]	57	SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )
[3]	58	!!----------------------------------------------------------------------
	59	!! * ROUTINE tra_zdf_imp *
	60	!!
[2602]	61	!! ** Purpose : Compute the after tracer through a implicit computation
	62	!! of the vertical tracer diffusion (including the vertical component
	63	!! of lateral mixing (only for 2nd order operator, for fourth order
	64	!! it is already computed and add to the general trend in traldf)
[3]	65	!!
[2602]	66	!! ** Method : The vertical diffusion of the tracer t is given by:
	67	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
[457]	68	!! It is computed using a backward time scheme (t=ta).
[2602]	69	!! If lk_zdfddm=T, use avs for salinity or for passive tracers
[3]	70	!! Surface and bottom boundary conditions: no diffusive flux on
	71	!! both tracers (bottom, applied through the masked field avt).
[2602]	72	!! If iso-neutral mixing, add to avt the contribution due to lateral mixing.
[3]	73	!!
[2602]	74	!! ** Action : - pta becomes the after tracer
[3]	75	!!---------------------------------------------------------------------
[2528]	76	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	77	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	78	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	79	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	80	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	81	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	82	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	83	!
	84	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	85	REAL(wp) :: zrhs, ze3tb, ze3tn, ze3ta ! local scalars
[6043]	86	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt, zwd, zws
[3]	87	!!---------------------------------------------------------------------
[3294]	88	!
	89	IF( nn_timing == 1 ) CALL timing_start('tra_zdf_imp')
	90	!
[6043]	91	CALL wrk_alloc( jpi,jpj,jpk, zwi, zwt, zwd, zws )
[3294]	92	!
	93	IF( kt == kit000 ) THEN
[457]	94	IF(lwp)WRITE(numout,*)
[2528]	95	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
[457]	96	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
[2602]	97	!
	98	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
	99	ELSE ; r_vvl = 0._wp
	100	ENDIF
[457]	101	ENDIF
[2528]	102	!
[2602]	103	! ! ============= !
	104	DO jn = 1, kjpt ! tracer loop !
	105	! ! ============= !
[2528]	106	!
	107	! Matrix construction
[2602]	108	! --------------------
	109	! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
	110	!
[2715]	111	IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. &
[2602]	112	& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
	113	!
	114	! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
	115	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt (:,:,2:jpk)
	116	ELSE ; zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
[2528]	117	ENDIF
[6043]	118	zwt(:,:,1) = 0._wp
	119	!
	120	IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
	121	IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
	122	DO jk = 2, jpkm1
	123	DO jj = 2, jpjm1
	124	DO ji = fs_2, fs_jpim1 ! vector opt.
	125	zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
	126	END DO
[2528]	127	END DO
	128	END DO
[6043]	129	ELSE ! standard or triad iso-neutral operator
	130	DO jk = 2, jpkm1
	131	DO jj = 2, jpjm1
	132	DO ji = fs_2, fs_jpim1 ! vector opt.
	133	zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
	134	END DO
[2528]	135	END DO
	136	END DO
[6043]	137	ENDIF
[2528]	138	ENDIF
[6043]	139	!
[2602]	140	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
[2528]	141	DO jk = 1, jpkm1
	142	DO jj = 2, jpjm1
	143	DO ji = fs_2, fs_jpim1 ! vector opt.
[2602]	144	ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	145	ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
[2528]	146	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	147	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	148	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	149	END DO
	150	END DO
[3]	151	END DO
[2528]	152	!
[2602]	153	!! Matrix inversion from the first level
	154	!!----------------------------------------------------------------------
	155	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	156	!
	157	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	158	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	159	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	160	! ( ... )( ... ) ( ... )
	161	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	162	!
	163	! m is decomposed in the product of an upper and lower triangular matrix.
	164	! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
	165	! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
	166	! and "superior" (above diagonal) components of the tridiagonal system.
	167	! The solution will be in the 4d array pta.
	168	! The 3d array zwt is used as a work space array.
	169	! En route to the solution pta is used a to evaluate the rhs and then
	170	! used as a work space array: its value is modified.
	171	!
[2528]	172	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
[2602]	173	! done once for all passive tracers (so included in the IF instruction)
[2528]	174	DO jj = 2, jpjm1
	175	DO ji = fs_2, fs_jpim1
[5600]	176	zwt(ji,jj,1) = zwd(ji,jj,1)
	177	END DO
	178	END DO
	179	DO jk = 2, jpkm1
	180	DO jj = 2, jpjm1
	181	DO ji = fs_2, fs_jpim1
[5034]	182	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
[2528]	183	END DO
	184	END DO
	185	END DO
	186	!
	187	END IF
[2602]	188	!
[2528]	189	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
[457]	190	DO jj = 2, jpjm1
	191	DO ji = fs_2, fs_jpim1
[5600]	192	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
	193	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
[6043]	194	pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
[5600]	195	END DO
	196	END DO
	197	DO jk = 2, jpkm1
	198	DO jj = 2, jpjm1
	199	DO ji = fs_2, fs_jpim1
[2602]	200	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk)
	201	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk)
[2528]	202	zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side
	203	pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
	204	END DO
[3]	205	END DO
	206	END DO
[457]	207
[2602]	208	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
[457]	209	DO jj = 2, jpjm1
	210	DO ji = fs_2, fs_jpim1
[2528]	211	pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[5600]	212	END DO
	213	END DO
	214	DO jk = jpk-2, 1, -1
	215	DO jj = 2, jpjm1
	216	DO ji = fs_2, fs_jpim1
[2528]	217	pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
	218	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	219	END DO
[457]	220	END DO
	221	END DO
[2602]	222	! ! ================= !
	223	END DO ! end tracer loop !
	224	! ! ================= !
[1438]	225	!
[6043]	226	CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwt, zwd, zws )
[2715]	227	!
[3294]	228	IF( nn_timing == 1 ) CALL timing_stop('tra_zdf_imp')
	229	!
[3]	230	END SUBROUTINE tra_zdf_imp
	231
	232	!!==============================================================================
	233	END MODULE trazdf_imp

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

Download in other formats: